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THE BROWN BOVERI 


The present article describes why Brown Boveri found it 
necessary to install computers, intended solely for technical 
and scientific services, and deals with the organization of the 
Centre in which they are housed. In addition to the a.c. net- 
work analyser there is an analogue and a digital computer. 
The article concludes with an account of some of the prob- 
lems of organization associated with the operation of the 


digital computer. 


The Need 


iy belong to the characteristic trends of modern 
technical development that the need of the 
engineer for means of simplifying calculation work 
should have grown at an almost unbelievable rate. 
There are various reasons for this phenomenon. One 
of the most important is the enormous increase in the 
unit capacities of the machines which have to be 
designed. In order to make the best use of the 
materials on the one hand, while avoiding the risk of 
difficulties later when the machine is in service, it is 
necessary to determine certain properties of the 
machine—for instance the forces produced by a 
short circuit, or the temperature rises in the copper 
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and iron—much more exactly than for a smaller 
machine. Moreover, in many cases the Company is 
obliged to give quite extensive guarantees, for 
example regarding the behaviour of a regulator in 
the event of a particular kind of disturbance. The 
ability to uphold such guarantees is often dependent 
on the results of preliminary calculation or experi- 
ments. 

Apart from the experiments during development, 
which continue to retain their full significance, the 
aids for calculation become increasingly important. 
The necessity for their employment primarily arises 
from the efforts made to improve the quality of tech- 
nical products. Of course cases may well be visuali- 
zed, particularly during the calculation work for new 
designs, where the ability to undertake the complete 
calculation of several variants allows the manufac- 
turing costs to be reduced. But, from the experience 
gained so far, this factor alone does not justify the 
purchase of a large computer, from the economic 
aspect. The development referred to above is taken 
into account by Brown Boveri. Since December 1959 


work has been in progress in a true Computer Centre. 
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Fig. 1. — Plan view of the Brown Boveri Computer Centre 


| = Network analyser room 
2 = Analogue computer room 


3 = Digital computer room 
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Fig. 2. — Analogue computer 


The various components are housed in the cabinets on the left. 

The switchboard on which the program for the particular cal- 

culation is set up can be clearly seen. To its right are an oscillo- 
graph and, behind it, a type KC voltage regulator. 


4 = Preparation room 
5-7 = Offices 


8 = Conference room 


This possesses three of the most important “‘calcu- 
lating machines’, namely an a.c. network analyser, 
an analogue computer and, as most recent but at the 
same time the largest member of the family, a digital 
computer which may be regarded as a calculating 
machine in the true sense of the word. The a.c. net- 
work analyser, which has rendered invaluable service 
for many years, is in effect a model of the network on 
which actual phenomena can be simulated. The 
analogue computer is also more of a model than a 
true computer. Admittedly the relationship between 
image and reality is less close in this case. 

Of course there are other mathematical aids in the 
factory, outside the Computer Centre, which may 
also be regarded as calculating machines, such as the 
electrolytic tank and the analogue model used for 
investigating the critical speeds of high-speed shafts 
Eye 4 


Layout of the Accommodation 


Since it was commissioned in 1953, the a.c. net- 
work analyser has been housed in the cellar of the 
main administrative building of the Company. It 
was therefore an obvious solution to locate the ana- 
logue and digital computers in adjacent rooms, there- 
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Fig. 3. — A.C. network analyser 


Panel with active and passive ele- 
ments and switchboard. On the left 
is the metering desk with automatic 


devices for calling up the elements, BROWN/BOVERI 


and light-spot instruments. 


by creating a Computer Centre. Fig. | shows the 
plan view of the Centre. The analogue computer and 
the network analyser (Fig. 2 and 3, respectively) are 
housed in adjoining rooms. With this arrangement 
it is possible to establish a connection between the 
two, as it is to be expected that problems will occur, 
demanding the combined efforts of the two machines. 
The digital computer (Fig. 4) also has a room of its 
own, in which the operator in charge of the machine 
can work undisturbed; but often the person con- 
cerned with the calculation being handled is also 
present. The punched tapes for programs and other 


Fig. 4. — Digital computer 


The computing and control sections 
with the magnetic-core memory and 
magnetic drum store are in the cabi- 
nets in the rear. On the left the 
punched-tape input equipment; in 
the middle is the operator’s desk with 
the output equipment (punching 
machine and page printer) on the 
right. 
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data are prepared in the adjoining preparation 
room. 

Naturally the operation of the Computer Centre 
needs a sufficiently large staff of experts, for whom 
offices are provided close at hand. Since the Centre 
has to serve a wide circle of “‘customers’’, a confer- 
ence room was also included. In it small groups 
can hold discussions, but it can also be used for 
internal lectures or courses of instruction. All rooms 
are air-conditioned, so that good and uniform 
working conditions are provided for the electronic 


computers. 
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Technical Data of the Digital Computer 


Reports on the analogue computer and network 
analyser have appeared in earlier numbers of this 
Review [3,4], so that further reference to their de- 
tails is superfluous. On the other hand, the digital 
computer (Siemens & Halske type 2002) deserves 
special mention. By modern standards this is a 
medium-sized machine. The model installed in the 
Company’s premises is ideal for technical and scien- 
tific duties. This means, among other things, that 
the computer itself, to judge by other electronic 
equipment, operates very rapidly, while the input 
and output sections are not quite so luxurious. Fur- 
thermore, for the applications involved it is impor- 
tant for numbers processed by the floating-point 
system not to be handled any slower or with any 
additional complications than they would by the 
fixed-point system, which is commonly used for 
commercial applications. 

Punched tapes are used for the input and output 
of information, facilities also being provided for the 
results of a calculation to be printed out in clear by 
a page printer. Advantage is seldom taken of this 
facility owing to the relatively slow printing speed 
of the page printer. For storing the program and 
data there is a magnetic-core memory capable of 
holding 2000 words and a magnetic drum for 10 000 
words. 

As regards the working speed of the computer it is 
not easy to make any binding statement, since this 
varies greatly from one problem to another. On the 
average, about 2000 operations (arithmetical and 


logical operations) are performed per second. 


Questions of Organization 


In conclusion it is worth saying something about 
the questions arising with respect to the organization 
of a Computing Centre of this nature to serve a large 
industrial undertaking. The field of application of 
a digital computer in the technical and scientific 
sphere in a company having the structure of Brown 
Boveri is extraordinarily broad, as may be gathered 
from the different subjects dealt with in the articles 
in this issue. Not only typical mathematical prob- 
Jems are dealt with (e.g. the solution of differential 
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equations) ; equally important are tasks which may 
not involve so much mathematical ability, but where 
series of routine calculations formerly carried out 
laboriously by hand can be handed over to the 
digital computer. 

There are two fundamental ways of running a 
computer centre of this kind; one is the “closed 
shop’’, the other the “‘open shop”. When the former 
principle is adopted, the problems to be handled 
by the computer are given to the staff of the Com- 
puter Centre in a concretely formulated form by the 
engineering departments. The Centre staff then carry 
out all tasks involved, until a numerical result is 
obtained. These tasks are: the determination of the 
mathematical method, the preparation of the block 
and flow diagram, encoding the instructions, check- 
ing the program and, finally, performing the calcu- 
lation according to the tested program, using the 
input data provided. 

In contrast, with the “open shop” system the 
greater part of the programming work is undertaken 
by the particular engineering department concerned. 
The Computer Centre then provides a technical 
service, acting in more of an advisory capacity; it 
enables the various department to exchange their 
experience, is responsible for the operation of the 
computer and—as one of its most important duties— 
ensures that its library of generally applicable sub- 
routines is utilized to the best advantage. 


Both systems naturally have their advantages and 
disadvantages. With the “open shop” system it is 
much quicker to spread the knowledge regarding the 
possible methods of electronic data processing in the 
engineering departments. Moreover, experience has 
shown that, for many tasks, the goal is attained much 
more quickly if an engineer of the specialist depart- 
ment is given instruction in the art of programming, 
than when a trained programmer has to delve into the 
special problems associated with a particular task. 
These advantages are offset by the undeniable draw- 
back that the programs produced by the open shop 
method are, on the average, less polished than those 
obtainable with skilled programmers operating as a 
closed shop. Weighing up the pros and cons it was 
decided from the very start to aim at operating the 
Computer Centre as an open shop. Several months 
before the computer was installed, internal program- 


| 
. 


MAY/JUNE 1961 


ming courses were held, so that now every engineer- 
ing department has at least one trained programmer. 
Experience gained in the meantime strengthens our 
conviction that this approach was right. Shortly after 
the “‘running-in period” had expired the computer 
was fully occupied; but this, of course, is no proof that 
the time spent in calculation was utilized to the best 
advantage. In the early stages allowance must be 
made for time being lost in checking the programs, 
owing to the lack of experience of the programmers. 
These are typical teething troubles, though, which 
can soon be overcome by well directed and rigorously 
executed measures. One of these, for instance, is to 
insist on the computer being handled only by the 
operator and not by the members of the engineering 
departments, even when checking through the 
program. Finally, experience has shown that a sui- 
table form of organization can always be found, even 
for special cases, provided development is followed 
with care and suitable precautions are taken. In the 


Brown Boveri Computer Centre the organization is 
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nearer the “open shop” than the “closed shop’’, 
although it exhibits certain features of the latter. This 
applies above all to scientific problems necessitating 
a very large amount of mathematical work. In such 
cases the Centre acts rather like a research depart- 
ment. 

M. CHRISTOFFEL 


(KME) 
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ON THE IMPORTANCE OF MATHEMATICS TO ENGINEERING 


With the aid of some typical examples, the present article 
discusses the mutual relationships between engineering and 
mathematics. Following a review of the pioneering days in 
electrical engineering, an attempt is made to find a means 


of determining the present state of the art. 


HE MATHEMATICIAN is, or at least until 

quite recently was, looked upon as an unworldly 
idealist, paying little heed to the possible applications 
for his scientific knowledge. This widespread im- 
pression is based on false assumptions though; for 
it is the great mathematicians in history who have 
repeatedly striven to place their ability at the service 
of physics or engineering. To mention but a few 
great names, readers are reminded of such men as 
Newton, Euler and Gauss. It is well known that 
Newton and Leibnitz were almost coincident in the 
creation of the infinitesimal calculus. By adopting 
new definitions Newton succeeded in formulating 
his famous laws of motion, which today are still the 
recognized foundation of the classical mechanics. 
Euler dealt with quite a number of technical prob- 
lems; he, for instance, is the father of the basic tur- 
bine equations. Finally Gauss—‘‘a prince among 
mathematicians’”—among his many achievements, 
rendered an invaluable contribution to the theory of 
electromagnetism. 

The number of mathematicians employed in indus- 
trial untertakings nowadays would have been beyond 
all imagination ten or twenty years ago. At the same 
time too, many engineers have experienced the need 
for expressing more and more technical relationships 
by mathematical means, with the accompanying 
trend towards a mathematical way of thinking. 
The fact that advantage can only be taken of the 
facilities offered by modern computers when there 


is a staff of skilled mathematicians available to 
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develop expedient numerical methods, is only a 
partial explanation for the above trend. We, rather, 
are more inclined to regard this as an. expression of 
the fact that modern engineering exhibits a pro- 
nounced trend towards the abstract. It is therefore 
most attractive to endeavour to trace this pheno- 
menon and, after a brief review of past history, to 
try to determine the present state of the art. In doing 
so, remarks will be principally confined to the field 


of heavy electrical engineering. 


Mathematical Thought 
in the Everyday Life of the Engineer 


When, in the remarks that follow, reference is- 
made to mathematicians and engineers, these two 
professions must be understood as types; it may well 
be that, in reality, an engineer possesses many of the 
qualities of a mathematician, and vice versa. But an 
engineer and a mathematician think on entirely 
different lines from one another. Those of the mathe- 
matician are centred on logic and the abstract, 
whereas the primary aim of the engineer is to create” 
something useful. Characteristic of this difference’ 
is the mental approach to mathematical formulae 
expressing the relationship between various physical | 
quantities. For the engineer a formula of this kind is’ 
generally little more than a convenient means of 
obtaining an unknown quantity from several known: 
quantities; the form in which it is expressed is. 
of little interest. The mathematician, on the other. 
hand, places great value on the aesthetics of a for= 
mula, not merely for the sake of the aesthetics as 


such, but because he primarily views the formula as a 


relationship between various quantities in accordance 
with a certain law, which he endeavours to express 
in its neatest and simplest form. 
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Now it is interesting to discover that, in the every- 
day life of the engineer, there are many indispensable 
conceptions which were created by mathematicians, 
but have become an integral part of the everyday 
language of engineers. Among these are, for instance, 


the terms function, variable, parameter, and so on, 


One of the great achievements of the philosopher 
Descartes was the discovery of analytical geometry, 
by which a link was forged between the geometrical 
conception and algebra. The very idea of expressing 
geometrical figures by algebraic formulae, on the 
one hand, and interpreting such formulae by geo- 
metry, on the other, was to have far-reaching conse- 
quences. Nowadays it is a matter of course for every 
engineer who has to perform a series of experiments 
to record the results, not merely as a table of figures, 
but also in the form of a graph, for which he utilizes 
the form of representation from which he can best 
obtain an empirical relationship between two or more 
measured quantities. In doing so he will begin by 
assuming that the desired quantity, termed mathe- 
matically as the dependent variable, varies in an 
unambiguous manner with respect to one or more 
other quantities, known as the independent variables. 
He will then perform a series of measurements in such 
a manner that for a particular group of readings only 
one dependent variable is allowed to change while 
all others are kept constant. He will attempt to 
plot the various readings belonging to one group in 
the form of a curve. Of course it will have to be 
assumed that the relationship in question is continu- 
ous, i.e. that a small change in one quantity does 
not lead to an abrupt change in the other. Often it 
is possible to pick out one of the dependent variables 
and to plot a complete family of curves on a single 
sheet of paper, on which the other variables now 
appear as parameters. He has thus succeeded in 
plotting the results of a long series of measurements 
on one graph which can be studied at a glance. These 
few remarks already give some idea of the tremen- 
dous amount of abstract thought which had to be 
carried out by the mathematicians of past centuries 


to enable this procedure to appear so self-evident 


days. 


Many other examples could be added, such as the 
e of complex numbers to represent alternating- 
nt quantities for calculations in electrical 
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engineering. Here we are not concerned with enumer- 
ating all possible examples, but rather with showing 
that even in technical spheres which appear to be 
extremely remote from mathematics, a certain 
amount of mathematical thought is involved some- 


where. 


Maxwell’s Theory of Electromagnetism and 
the Early Days of Electrical Engineering 


In the first half of the 19th century countless 
attempts were made by the mathematicians and 
physicists to summarize the numerous discoveries 
regarding electromagnetic phenomena by means of a 
self-contained mathematical theory. Not one of these 
attempts led to a conclusive success, and they are 
now all forgotten. It was Maxwell who succeeded 
in postulating the theory which found its final form 
in his now famous “‘treatise on electricity and ma- 
gnetism”’ published in 1873. Aided by the experimen- 
tal work of the renowned Faraday, with his lucid 
methods of illustration, he managed to summarize 
all the electromagnetic phenomena in four simple 
equations. If these are augmented by the three “laws 
of matter”, i.e. providing the relationship between 
electrical field strength and dielectric displacement, 
or between magnetic field strength and flux density, 
together with Ohm’s law, it is possible to explain 
all the phenomena known at that time, with the aid 
of this theory. The postulation of this theory is justly 
hailed as one of the greatest scientific achievements 
of all time. It is characteristic of a successful theory 
that, like Maxwell’s, the equations not only satis- 
factorily explained electromagnetic phenomena 
known at the time they were postulated, but also 
permitted prognoses regarding phenomena which 
were still unknown. For example, Hertz discovered 
radio waves on the basis of predictions from Maxwell’s 
equations, and the experimental discovery that the 
so-called electromagnetic constant is identical with 
the speed of light, led to the realization of the fact 
that optical phenomena involve electromagnetic 
waves. 

A further feature of Maxwell’s theory is its des- 
criptiveness. The first two equations in their integral 
form permit of a concrete interpretation and are 


nowadays known to every high-school boy as the 
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law of induction and the law of magnetic flux 
(‘“Maxwell’s Rule’). The two other equations ex- 
press the fact that, firstly, the paths of electric current 
are always closed and, secondly, that the source of the 
electric field must be sought in the electric charge. 
The lucidity of this theory—consider for instance 
the terms inductance, capacitance, reluctance, leak- 
age factor—is certainly one of the main reasons for 
the victorious advance made by electrical engineering 
at the end of the last century and the beginning of 
the present one. To understand the operation of 
rotating electrical machines and transformers there 
is no need for any special mathematical knowledge, 
and they can also be designed with simple formulae. 
It is therefore obvious that, at that time, intuition 
and imagination were more useful to the engineer 
than special mathematical knowledge. Even so, there 
were certain problems which necessitated intensive 
development of Maxwell’s equations. One of the 
most notable is Rogowsky’s calculation of the leak- 
age field of a transformer. It is proof of Rogowsky’s 
great mathematical skill that he should have suc- 
ceeded in obtaining a final result, after considerable 
effort involving the expansion into Fourier series, 
leaving only one unexplained factor which is now 
known as the Rogowsky factor. Another phenomenon 
is the ‘‘skin effect’? which causes the additional 
losses in the windings of electric machines. In this 
case, Field, starting from Maxwell’s equations, postul- 
ated a very elegant theory, which allows these losses 
to be predetermined quite easily. 


Modern Trends 


Without exaggeration it may be claimed that 
mathematical knowledge and methods occupy a 
more important position in engineering than ever 
before. It is possible to give a number of reasons 
explaining this phenomenon, which will now be 
briefly outlined. 

Naturally, the most obvious is the fact that 
nowadays, with electronic computers, it is possible 
to tackle problems which would formerly have 
been a hopeless task for an engineer equipped only 
with slide-rule or, at the best, a portable calculating 
machine. But it is a well-known fact that the elec- 
tronic computers only execute the orders they are 


given; therefore suitable numerical methods are re- 


quired to enable a computer to solve a complicated 
differential equation, for instance. Numerical analy- 
sis, as a branch of applied mathematics, has conse- 
quently witnessed a tremendous boom and fills an 
important place in the curriculum of universities. 

It is remarkable that the need to make use of 
these new facilities is also felt in the typically 
classical spheres, such as heavy mechanical and 
electrical engineering. This is quite comprehensible 
if we consider the size of the machines involved 
nowadays. The enormous increase in the unit outputs 
of steam turbines, generators and transformers is not 
only due to technological progress; it would not 
have been feasible if the methods of calculation had 
not been improved at the same time. Only in this way 
has it been possible to achieve the double goal, of 
increasing the power-weight ratio while improving 
the reliability. It is quite obvious that, for instance, 
the mechanical problems associated with the cons- 
truction of a 300-MVA generator need studying in 
far more detail than those for a machine with an 
output of the order of 10 MVA. 

In the research programme of most engineering 
firms the study of transient phenomena of all kinds 
occupies an important position. In this field in 
particular the demands have increased very greatly. 
There was a time when it was regarded as a 
natural occurrence for a transformer to break down 
as the result of a lightning stroke. Nowadays it is 
taken for granted that a transformer will be so 
designed that it can withstand such stresses (provided, 
of course, it is equipped with a suitable means of 
overvoltage protection). Furthermore, it is remark- 
able that a multiple winding system, as represented 
by a transformer, possesses such a complex con- 
figuration that nobody has yet succeeded in estab- 
lishing a theory for transients which proves satis- 
factory in all conditions. 

Also in rotating machines transients, i.e. their 
behaviour when the load suddenly changes, are 
important. Naturally, generators were being built 
long before any theory existed regarding transients. 
Even today this theory has not been fully developed’ 
in all details. In this period of growing unit ratings 
and increasing interconnected operation, however, 
there is an obvious need for continuous refinement 
of the theory and methods of calculation. 
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The ability to perform extensive calculations 
quickly provides a notable incentive for theoretical 
work. Not so very long ago it was possible to speak 
of stagnation in this field. The basic operations were 
known in principle, and refinement of the theory 
was pointless, so long as there were not any numerical 
methods of utilizing it in practice. Compared with 
earlier eras, it is striking how short a time it now takes 
for new theoretical knowledge to become general 
knowledge among engineers. It may be recalled, for 
example, that Steinmetz and others recognized before 
the beginning of the present century that steady-state 
alternating-current quantities could be represented 
by complex numbers which could be employed for 
calculation. On the other hand, a text-book on 
electrical engineering! published in 1922, which had 
the reputation of being very up-to-date at that time, 
did not even mention this ability. Thus the present 
situation is stressed even more clearly. 

What has been said so far in this article regarding 
spheres which may be regarded as classical, applies 
to an even greater extent to more recent technical 
creations. Whereas it might, perhaps, be feasible 
in an emergency to continue building generators and 
transformers without any improvement of the basic 
theory, or increased calculation, this is certainly 
not true of such fields as nuclear physics or guided 
missiles. Here the outlay for experimental develop- 
ment is so great that technical development is quite 
out of the question without well-founded mathe- 
matical theory and the facilities for carrying out 
extensive numerical calculations within a reasonable 


space of time. 


1 Theoretisches und praktisches Lehrbuch fiir Elektro- 
techniker. J. Fischer-Hinnen, Zurich 1922. 


At this point it is certainly fitting to damp some 
of the enthusiasm regarding the capabilities of 
electronic computers. There are a great many prob- 
lems whose numerical solution involves such a large 
number of operations that even the fastest electronic 
computers cannot manage them. We are mainly 
considering the solution of partial differential equa- 
tions with more than two independent variables. 
Of practical importance among such problems is the 
calculation of three-dimensional fields where, even 
today, one cannot always be certain as to how these 
can be calculated within a reasonable space of time. 

Finally it is worth saying a word about a branch 
of mathematics which has made great progress in 
the last few decades, namely statistics. The applica- 
tions of this science continue to gain significance in 
engineering. The reason is obvious. Extended series 
of experiments, e.g. breakdown tests on special 
models, mostly take a long time and involve heavy 
costs. By means of statistical methods it is possible, 
with good advance planning, to gain a maximum of 


information with a minimum of experimental effort. 


Conclusions 


If we follow present technical developments 
attentively, we are bound to arrive at the conclusion 
that ‘‘mathematization”’ is making rapid progress in 
engineering. This implies that higher standards have 
to be set for the theoretical and mathematical train- 
ing of engineers. But of course it does not mean that 
they can dispense with the main prerequisite for 
the creation of new technical products—an inventive 
imagination. 

M. CurisTorFEL 


(KME) W. Frey 


292 


THE Brown Boveri REVIEW 


VOL. 48, No. 5/6 


THE INTEGRATION OF FUNCTIONS 
WITH THE AID OF ELECTRONIC COMPUTERS 


For the evaluation of simple and multiple integrals the 
trapezoidal rule is the simplest method. However, it is seldom 
used because it is prone to introduce a much larger error than 
other methods, such as Simpson’s formula. On the other hand, 
it offers the advantage that the calculation with half the 
integration step can be solved without difficulty, there being 
no need to maintain the functional values of the first calcula- 
tion for the double step, or for them to be calculated again. 
Extrapolation in terms of h? according to Richardson (where 
h denotes the length of the step) yields the same results as 
Simpson’s method when the trapezoidal rule is used twice 
in succession, with steps in the ratio of one to two. Thus the 
two methods are equivalent. By utilizing the Richardson 
extrapolation once more it is possible to determine whether 


the accuracy attained is adequate. 


HE INTEGRATION of functions is one of the 

fundamental operations of mathematics. For the 
solution of many problems numerical integrals be- 
tween wide limits are required, or multiple integrals, 
for example for the calculation of the forces produced 
by currents flowing in systems of conductors. 

‘To attempt to undertake all the calculations needed 
is far beyond the scope of human capabilities, so that 
electronic computers have to be employed. For this 
it is necessary to employ definite methods, which also 
have to be tailored to suit the particular case of 
utilization with electronic computers. The question 
of the accuracy of calculation is closely related to the 
choice of the integration interval. In the course of 
this article these relationships will be illustrated by 


reference to an example. 


Peculiarities of the Electronic Method 
of Calculation 
The technique and capacity of numerical cal- 


culation was radically changed by the introduction 
of desk calculating machines. This applies to a much 


681.14: 517.3 


greater extent to the application of automatic com- 
puters. 

An outstanding feature of the special methods used 
with electronic computers is their simplicity com- 
pared with some of the relatively complicated integral 
methods put forward by certain authors [1, 2, 3]. 

The following remarks are primarily intended to 
outline some of the advantages of the simpler 
methods, especially that using the trapezoidal rule: 

Leu~- { 
. b-—afl MEX 
Uide came Eee: +9) + y">,| 


4 t=] 


(1) 


with 3, = f(a +i" be 4 


Above all, this formula is easily programmed; 
there is no need to store any constants in the com- 
puter. Also the number of integration steps n— 
within the total interval from a to b—enters into the 
calculation in a very simple manner. Particularly 
when a wide integration interval is involved, the 
constancy of the step length exerts a_ beneficial 


influence on the work of programming. 


As will be known, the error is proportional to h2, — 


where h = (b — a)/n denotes the integration step. 
In the next chapter a method will be demonstrated 
in which, by means of an artifice, it is possible to 
attain the same convergency, proportional to h*, as 
is attainable with Simpson’s formula. 

The following reason is given for the preference 
shown to the trapezoidal rule. Assuming the integra- 


tion were carried out with a step of length h, where — 
vie is easy to keep in store, and the calculation — 


were subsequently repeated with half this step 
length, only the values in the middle of each step 
of the first computation need to be calculated, the 
values so obtained rendering a contribution to those 


. 
F 
q 


» 
b 
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determined with a step of half the length. This 
advantage is particularly marked when the calcula- 
tion of f is tedious, especially when double integrals 
written in the form 


j (‘i (xy) dy Jas 


a \p(x) 


have to be dealt with. This argument applies even 
more forcibly to multiple integrals. 


Extrapolation in hm 


In 1927 Richardson put forward a method of im- 
proving the accuracy of numerical results, provided 
the order of magnitude of the error is known [4]. 
The basic idea of this method is easy to under- 
stand. 


It may be assumed that, if A is sufficiently small, 
the error is proportional to h”. Furthermore, the 
calculation is considered to have been performed with 
two values of A, e.g. with h; = h and hy = 2h. If we 
desire to determine the exact value J, we must assume 
that the following values J; and J; are known in 
advance. 


fy 14K hi™ 


(2) 
Tete Rh” 


By eliminating K we obtain 


ht i ; (3) 

hy 
A prerequisite condition for this result is that the 
relationship obeys equation (2), which applies for 
small values of h. The application of formula (3) is 
limited by round-off errors. Nevertheless it can be 
checked by calculation whether the equation (2) 
is fulfilled, by calculating J with (3) using different 
steps and comparing the results. 

In the example given below the case is considered 
for n= 4 and n= 2, with hj = 2 hy; the latter 
assumption implies that each time the integration 
is performed with the step reduced by half, thereby 

facilitating the execution of the calculation. 


For an operation with an error proportional to h? 
we obtain 


4 OT, 
(=a —_— 3 
j (3') 


For an error of the order of magnitude of ht 


ib ieee 
oa TEE 


(3") 


At an earlier stage it was mentioned that the trape- 
zoidal rule yields an error of the order of h?, hence 
the equation (3’) is valid. 

The trapezoidal rule (1) is now applied to the 
integration with two different steps A and 2h. Then, 
with n = (b — a)/h it becomes 


b—a({l 1 | 
Ed se ie qt RtM A tI + oI 
r) 


b—a{l 
= Oy ay eee IR aie eae 


n 
] 
Siete + Jn-2 + Jn-1 =e 9 »| 
From equation (3’) we obtain 


aeee 


b 
1 Torrpres ee che BAN aie te ee eR Ye 
Ne aa 2 Yn-2 ah 4 nt + Jat (4) 


This form is identical with Simpson’s formula, for 
which the error is known to be of the order of At. 
Supposing S, is the result of the computation with 
the step A, then 


_ 4 iF a To, 


Ss 
s 3 


(4’) 
It is also possible for the values obtained from 
Simpson’s formula with the partial interval h to be 
extrapolated with f”, as follows: 


pa 16 Si, jee Soy, 


I 


From this aspect extrapolation in hk” automatically 
allows the length of the step to be controlled. Another 
possibility is afforded by the iteration of S,, S,, Sjy2 
and calculating J each time. If, in the course of this, 
the difference between two successive values of J is 
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less than a given amount ¢, the last value obtained 
is regarded as an approximation for the desired 


integral. 


The procedure is as follows: 


Starting with h = (b — a)/2, the following values 
are calculated with the aid of the trapezoidal rule, in 
that the number of steps in the integration is doubled 
each time. If ¢, denotes the value obtained with 2” 


intervals, we obtain for n > 4 


64 t, — 20 ti + the 
= " a 
n 45 ( ) 
This formula corresponds to (5) in which, with the 
aid of (4’), the value given by Simpson’s formula was 
expressed by those obtained with the trapezoidal 
rule. The calculation is discontinued when J, — J, is 


smaller than «. 


Example 


The starting point for the above study was the 
consideration that it would be illogical in the case 
of a wide integration interval for points situated 
within the limits to be differently weighted. Thus, 


with the previous notation, the following formula is 


This formula provides accurate results for poly- 
nomial expressions of the third order. Since it is 
extremely complicated to estimate the error, evalu- 


ation is replaced by an experimental method 
10 
Supposing, for example, i) sin x dx = 1-8390715 
0 


were to be calculated. 


In Table I the values obtained by the three 
methods given above—trapezoidal rule, Simpson’s 
formula and equation (6)—are tabulated to the 
seventh place of decimals, enabling them to be 
compared. 

It will be recognized at once that Simpson’s 
formula yields more accurate results, i.e. smaller 
error, than the other two methods. The order of 
magnitude of the error can be estimated as follows. 
Supposing J is the exact result and, with the aid of 
the trapezoidal rule, we obtain 


I=t,+ Kh? whereh = (b —a)/n 
and P= tae] +R 


Eliminating K we obtain 


[Se 
IT = tay 


For a formula of the order h” the relationship be- 


obtained : : : ete 
tween two successive differences, in the limiting case, 
: b—a ee is 2”. 
idxi— : : : 
J - n a +1 (Jo +In) + 2 ts Table II shows the ratios of successive errors in 
(6) terms of n. 
TABLE I 
Number of Formula (6) Simpson’s formula Trapezoidal rule 
partial 
intervals Value Error* Value Error* Value Error* 

2 —7-2995303 91386018 —7+2995303 91386018 —6:1546741 79937456 

4 0-9961065 8429650 3-0700157 12309442 0-7638433 10752282 

8 16615958 1774757 1-8695611 304896 1-5931317 2459398 

16 1-7964553 426162 1-8407061 16346 1-7788125 602590 

32 1-8285228 105487 1-8391702 987 1-8240808 149907 

64 1-8364410 26305 1-8390778 63 1-8353285 37430 

128 1-8384144 6571 1-8390721 6 1-8381362 9353 

256 1-8389074 1641 1-8390717 2 1-8388378 2337 

o12 1-8390306 409 1-8390717 2 1-8390132 583 

1024 1-8390614 101 1-8390717 2 1-8390570 145 


* Given as absolute values x 10? 
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TABLE II 


Trapezoidal Formula (6) Simpson’s 
rule formula 
7:4344642 10-841022 7-4240585 
4-371916 4-7497488 40-3726 
4-081349 4-16451 18-653 
4-01976 4-03994 16-56 
4-0050 4-0102 15-7 


4-0019 4-0032 _- 
4-002 4-012 


4-02 4-05 
TABLE III 

n Formula (5’) Simpson’s formula 
8 1-789531 1-869561 
16 1-828782 1-840706 
32 1-839068 1839170 
64 1839072 1-839078 
1-839072 1839072 


Exact result: 1-+839072 


From the above Table II it is evident that the 
quotients tend towards a limit, the round-off error 


becoming increasingly pronounced. For both the 
trapezoidal rule and the formula in equation (6) it 
attains the limit 4. This indicates that for both 
methods of calculation the error possesses the order 


of magnitude of h?. The theoretical proof of this fact, 
however, encounters considerable difficulties. 

In conclusion the trapezoidal rule is now compared 
to the expression (5’). Table III shows an improve- 
ment on Simpson’s formula, in that it manages with 
one iteration less. 


Conclusions 


The examples listed in the Tables I to III yield the 
following results. 

Formula (6) should not be used as it does not pro- 
vide any better results than those obtained with the 
trapezoidal rule. Experience confirms the correct 
estimation of the error, thus justifying the use of the 
suggested method of calculation. Hence useful results 


may be expected with this method. 


(KME) P. BANDERET 
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THE USE OF DIGITAL COMPUTERS TO CALCULATE BESSEL 
AND ALLIED FUNCTIONS 


The article gives the polynomial approximations for modi- 
fied Bessel and Struve functions and their integrals for the 
parameters 0 and | in the interval 0 < x <o, which are 
often required for the computation of axially symmetrical 
problems such as those connected with the magnetic fields 


on transformers with concentric windings. 


OR THE AUTOMATIC calculation of trans- 

cendent functions by electronic computers, the 
best expressions from the programming aspect are 
usually provided by cut-off power series or asymp- 
totic series. For large arguments, though, power 
series may be numerically unfavourable, owing to 
the round-off error. On the other hand, sometimes 
the asymptotic series cannot be used with the same 
x values (e.g. for Ko(x) ), because the accuracy there- 
by obtained is inadequate. Generally speaking, 
another method is more appropriate. In that case 
polynomials of the form P,(x) are employed which 
afford a close enough approximation to the function 
f(x) in a given interval (xo, x,). Thus for these 
polynomials 

Max | f(x) — P,,(x)| 


xy Sx Sry 

should be smaller than 
Max |f(x) — P,, (x)| 
4 SES, 


in which P,(«) is any polynomial whose order does 
not exceed n. The polynomial P,(«) is known as the 
Chebyshev approximation of the n-th order for 
f(x). A feature of it is that the maxima, of which there 


are not more than n + 2, for the expression 
| f(x) — P,(x) | 


are all equal in the interval (xo, x) [1, 2]. 


681.14:517.5 


Since the calculation of the coefficients of P,(x) is 
extremely tedious, a number of numerical methods 
have been devised, which do not yield the Chebyshev 
approximation itself, but provide polynomials of 
almost the same accuracy. Two of these methods, 
which have given satisfactory results, will be utilized 
in the following remarks. Other feasible methods, 
such as the one utilizing rational approximations, as 
propounded by Maehly [4, 5], are not considered. 

During the calculation of the magnetic field of 
transformers with concentric windings it became 
necessary to calculate Bessel and allied functions, for 
which great accuracy was stipulated, and the shortest 
possible time allowed for calculation. For this task 
auxiliary polynomials were determined; with their 
aid and the elementary functions e*, log, « and Vx, it 
is possible to calculate the modified Bessel functions 
Io, Ii, Ko, Ki, Jo —Lo, 1; -—L, (£ being the modified 
Struve functions), as well as the integrals fx I, dx, 
fx K, dx, fx(Z, —Ly) dx. The accuracy of the 
the entire interval 


auxiliary polynomials over 


0 <x < o, amounts to nine significant digits. 


Method of Calculation 
The polynomial approximations for 
Ko, Ki, [x Ky dx 


were determined, using the Lanczos t method [1, 3]. 
Attempts to adopt the same method for 


Ih, hh, [x Ly dx, Ip—Lo, h—-Ly, [x(hy — Ly de 


i 


proved less successful, so that, for these functions, 
approximations according to the telescoping method 
were calculated. 


a 
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The t Method 


It is assumed that the function » = f(x) is given 
as the solution of a linear differential equation 
L(») = 0. The coefficients are assumed to be poly- 
nomials in x. Lanczos’s t method now permits a con- 
vergent series of polynomials of increasing order to 
be formed, in which the accuracy improves with the 
order of the approximation polynomial. 

For this purpose the equation L(y) = 0 is con- 
verted so that a term TZ", (x) has to be added on the 
right-hand side. tT is a small, constant number and 


/. 


(x) is the Chebyshev polynomial of the n-th 
order in the interval (0,1). The resultant, modified 
equation 


L(y) =L (9) — tT, (x) =0 


then yields a polynomial as solution, and this exact 
solution of the modified differential equation is an 
approximation for f (x) in the interval (0,1) with 
practically the same accuracy as the Chebyshev 
approximation polynomial of the same order. 

The t method often also yields approximations in 
the interval (a, cc) in those cases where the asymp- 
totic series is inapplicable; on the other hand, it is 
not always possible to predict whether the method 
will yield usable results or not. 

An example will now be shown in which the 
method is employed twice in succession. In order to 
obtain a series of convergent polynomials in the 
interval 1< x < ~ for the function 


r/2 — [Ek (&) dé 


it is necessary to start from the asymptotic form of 
the function K,(x) 


K,(x) =Felet St -)= TE 2(3) 


The function z(1/x) then satisfies the differential 
equation 


ey) op im 


1 The constant is chosen so that the function tends towards 
0 when x+ oo [8]. 


obtained from the differential equation for K,(x) 
org l 
when Vz Zz (4) is substituted therein for K,(x). 


With the formal development 


and adding t, T, (1/x) to the right-hand side, the 
following recurrence formula is obtained for the 
values of b, (k =n, ... 0) 


— Fille — 8 (kK + 1) bet (3) 
: (2k + 1)2 — 4 


in which ¢,,, 1s the coefficient of the k-th power of 1/x 
in the polynomial T,, (1/x). 


With the supplementary condition 
bn41 es 0 


which follows from (2), and the further condition 


by = Vr/2 
resulting from the boundary condition for K(x) at 
infinity, the coefficients 6, can be unambiguously cal- 
culated from (3). The above supplementary con- 
dition 5, ,; = 0 also fixes t. This method can be used 
when t becomes sufficiently small for a definite 


value of n. 


The most obvious asymptotic series which can be 
used to represent the function 


x 


z/2 — | & Ky (6) dé 
is obtained by partially integrating the asymptotic 
expansion of the function x K,(x), term by term, 


taking into account the boundary conditions 


j ~« Vx o = 
n/2 — JE Ky (6) dé=e V (+ 4 +...)= 


——J e* x: y fe 


The function y (1/x) then satisfies the differential 


equation 


Fb) 


when z (1/x) is the function defined by (1). 
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With an expression for y,_; (1/x), analogous to that 


for z, (1/x) in equation (2) 


] a\ 
Yet{—}= a9 +—+.-.+— 
xy x 


a further recurrence formula is obtained for the 
coefficients a, (where k=n-—1,..., 0) from (4), a 
term T2 JT, (1/x) having also been added to the right- 
hand side. This formula is 


T2 bnkt1 — 2 aey1 — 2 beg (5) 
2k—1 


Since the conditions 


a,= 0 and a = 3 


apply here too, and the values of b, are known from 
(3), it is possible to calculate the coefficients a, 


unambiguously. 


Telescoping Method 


Supposing a power series 
P(x) = a9 + ax pe aX (6) 


is given, which represents the function y = f (x) 
sufficiently accurately in the interval (0, 1), the last 
coefficient a, should not be below the limit of the 
accuracy, i.e. the term a,x” is not so small that it can 
be neglected. The telescoping of P, (x) is now based 
on the approximation of all powers of x by polyno- 
mials of lower order, but with practically the same 
accuracy. The coefficient ¢,, in the Chebyshev 
polynomial 


T,(*) = ban X” + tant Rez? +... tn,0 


is equal to 24”, so that 


es tee ee 
O21 = ae es 
or, also, 
Te 
= — (apy xe! +... +.) 4 ae 


Thus if we substitute 


— (%-, x71 + ... + a) 


for x” in (6), we introduce an error which, in the 
entire interval (0,1), at no time exceeds the value 
of 2 x 4", and the calculation of P,(x) is reduced 
to that of a polynomial P,; (x) of the (n — 1)-th 
order. This process of telescoping can now be repeated 


until the last coefficient a,, is so large that even 
multiplying it by 2 « 4” cannot diminish it enough. 
The telescoping process is then finished and the 
original series, an approximation to the Chebyshev 


approximation, determined. 


Approximation Polynomials for 
the Modified Bessel and Struve Functions 
and their Integrals 


The approximations given below have been cal- 
culated in such a way that, for the interval (0, a), the 
corresponding power series were improved by the 
telescoping method. For the interval (a, co) either 
the corresponding asymptotic series were truncated, 
also by the telescoping method, or the t method was 
utilized to obtain the approximations. In the process, 
attention was paid to preventing the error from ex- 
ceeding five units of the tenth significant digit over 
the whole interval (0, ). 


il Io (x) = Io 
x \2k 
ORerrcanl\() i lp € 
0 SS ( j a oa 
k 
e* 10 
10 <tc cope eee ye (-) dire 
Vx 1% Z 
ay | 
Spe 0-00000 0003 0-00000 00002 
k=0 1-00000 00027 0:39894 22804 7 
1 24-99999 90983 0-:00498 67722 4 
2 156-25005 02179 0:00028 06006 0 
3 434-02668 00860 0:00002 86977 0 
4 678:18085 14731 0-:00000 58039 7 
5 678-08454 24100 —0-00000 07275 3 
6 471:31196 73418 0:00000 10787 1 
i 239-23918 92805 
8 95-90195 23737 | 
9 26°24210 75769 
10 9-66646 25536 
11 0:14399 45401 
12 0-66883 15150 
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y = 0-57721 56649 01533 


> 


0,00000 00001 


0 0-50000 00000 9 
1 6°24999 99695 0 
2 2604166 83656 6 
3 54+25343 50705 6 
4 67-81726 17980 4 
5 56°51119 24416 | 
6 3365156 77364 2 
7 14-98220 50575 7 
8 5-28380 36833 2 


9 1+35566 93198 1 
10 0°41187 05935 5 
11 0-01604 45909 7 
12 0-02411 17429 2 


0,00000,00002 


0-39894 22803 2 
—0:01496 03283 8 
—0:00046 76177 3 
—0-00004 03143 9 
—0:00000 72622 7 

0-00000 07463 2 
—0:00000 12354 5 


3. i EI, (6) dé = i, 


Drax <= 19 ia |Ya 


x 


0-00000 00000 5 


0° 16666 66666 97 
4-51249 99781 47 
48-48251 74891 78 
283-60014 30398 28 
1047-06815 64989 18 
2665-27067 48766 90 
4964-09851 71536 51 
7054-95843 37915 41 
7933-29573 08716 29 
7116-15292 07114 55 
5571-34195 40440 09 
2960-52407 13922 96 


o 
| 3 
7 
* 


OoOmnroauhowond— © 


12 2532-12381 56837 46 
13 ~94-65220 30647 81 
14 1245-53338 57226 89 
15 —533°52266 96117 54 
16 366-76444 46018 50 
17 —97-36153 67903 26 


21-24575 56155 02 


4. Ko(x) = Ko 


ita l 4 
= (2 OB 5): 


0-00000 00003 


0-39894 22804 2 
—0-01837 23441 5 
—0-00061 29416 9 
—0-00005 46931 2 
—0-00000 64035 4 
—0-00000 36450 6 

0-00000 17385 3 
—0-00000 09470 5 


5 ek + b el 
de | y, k | 


b xix no! Ks =| da(z) +4| 


af 


0-00000 00007 


—1-00000 00001 3 
1 | -0-24999 99933 8 
—0-01562 50529 8 
—0-00043 38794 3 
—0:00000 69512 3 


Ke COMO WAND Of OF ND 


— 


k 


by 


0-00000 00000 1 


0-25000 00000 02 
002343 74998 87 
0-:00479 57184 94 
0-00001 41259 84 
0-00000 01577 33 


Ck 


0-00000 0002 


1+25331 41373 
—0-15666 40293 
0-08809 47319 
—0-09099 77576 
0-13035 95356 
—0-20650 10623 
029703 31052 
—0-33905 07775 
0-28091 80905 
—0-15577 95499 
0-05128 12921 
—0-00754 46371 


Sere a 


1 
(ikea cen | R= t(7 + lows 


k 


Ky 


2 


[Nawttel+s 


x 


b, x7* + ¢ 
het 


y = 057721 56649 01533 


e- 
|| ey eS eve hey = — 


se 


a 


Emax | 0*00000 00000 5 


k=0 |0-50000 00000 22 
1 |0-06249 99988 96 
2 |0-00260 41754 97 
3 |0-00000 42287 47 
4 |0-00000 07064 25 


k 


ri 


0-00000 00000 9 


—0°25 
—0:07812 5 
—0:00434 02777 


—0-00010 62463 83 


—0-00000 14806 68 
—0-00000 00133 75 


| 


ok 


0-00000 0002 


125331 41373 
0:46999 25204 
—0-14683 86549 
0-12757 69501 
—0-16879 46725 
025702 30979 
—0-36241 23971 
0-40933 35460 
—0-33710 41487 
0-18623 64530 
—0:06115 49277 
0-00898 15841 
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6. x/2 — [EK (€) dé =hy 


O<xe< 1 he=a/2—-* 7% (y+ lo vl 


0 


a, 


0-00000 00000 4 


: k=7 80546-90721 62114 48 

8 1 91082-95533 82325 02 

| Soa a2* + | + x? paz xk | 9 | 4 03861-79309 34627 77 

i v 10 7 64503-93519 33666 21 

11 | —12 95216-38659 23775 26 

y = 057721 56649 01533 12 | 19 50764-97690 44749 25 

13 | -25 81260-86233 18439 89 

= 1\* 14 | 29 54256-88282 76632 76 

l<x<co ky=e*/x [Dal =) ae | 15 | -28 71496-73890 00573 20 

: 16 | 23 21565-07609 97806 52 

17 | -15 23703-71975 31415 93 

ie e : 18 7 87453-31186 20749 97 

19 | 3 07281-58366 25414 90 

20 84863-78832 02513 07 

Emax | 000000 00000 05 [0-00000 00000 08 | —0-00000 0002 as SPE My eet 

k =0| -0-16666 66666 675|0-13888 88888 910} 1-25331 41373 22 1211-58770 94691 98 
-0-01249 99999 575|0-01812 49998 963) 109664 88851 
-0-00037 20241 492|0-00067 31867 710] -0-69509 10682 


—0-00000 60272 126 


0-00001 24726 259 


1-16661 91735 
—2-99524 18887 
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0-00000 00007 


0 

2 

3 

4 |—0-00000 00627 385\0-00000 01428 666 
5 9-43322 66362 
6 

7 

8 


8. 1, (x) — £, (x) = S 
3061340 88398 TA aes : 
9041176 79352 : 
“4 
22605985 97708 0.< x29 Se Y" a, =) eee 
9 460-46962 68423 , 19 
10 -747-44666 97296 
11 951-10275 99993 2 ] 19 \2 
19< 2 ee Feset 
12 9329-46524 38657 #5 CS nly Ui aes | b ( pee 
13 688-19980 22368 
14 -369-07796 01792 
15 135:57552 40815 
16 -30-47957 22146 a | by 
17 3-16059 58754 
Sak 0-00000 00000 3 0-00000 03 
i= 0:49999 99999 79 0:63661 9700 
1 4.03192 52059 27 000529 2639 


2 2256249 71324 95 
3 —97-03482 85493 96 
4 339-37212 81780 86 
5 1000-73942 09260 89 
6 
7 
8 


0-00020 9945 
0-00003 4520 


fx» k 
0 <x < 20-25 i Ne a( lee 
d 20-25 
1 20-25 \ 2h 2550-99133 94458 52 
20:25 < x < 00 So=y Yi b( ‘3 zig © 5721-81554 06125 88 
s 11426-34647 35352 31 
9 | ~20429-85677 72348 26 


a | - 10 32659-60240 77395 51 

11 | 46278-37659 26248 03 

12 57253-74660 14571 80 

Emat 0-00000 00000 4 0-00000 00007 13 60622-87018 55409 99 
k=0 0-99999 99999 65 0-63661 97726 9 14 53665-13699 96938 17 
1 ~12-89155 03526 05 0-00155 24793 6 15 | -38656 42330 24412 35 

2 102-51561 81809 52 0-00003 41867 7 16 21927-84744 97397 40 

3 ~587-37077 45667 79 0-00000 17952 6 17 -9377-30090 93586 30 

4 2627-34459 14092 55 0-00000 05069 7 18 2830-60011 95461 30 

5 ~9633-91509 95177 90 19 —536°00645 34299 52 

6 29920-29683 84453 68 20 47-78333 21376 74 


ee ee | 
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9. fé [1 (€) — L1(6)] dé = 5 


k 
| a Wy Ae =#|Ya(77) + | 
re 


<a xX< co 5S) = ¢ + 


NG 
1 17 \24 
pd ee € 
¢ = — 0-17211 83100 
ag | by 
Emax 0-00000 00000 1 0-00005 

e— 0 0-16666 66666 59 0-95495 205 

1 —0-90187 80049 57 002441 489 

2 3-61249 92739 70 0-00282 340 
3 —11-58408 64250 65 
+ 31-07087 55157 04 
5 —71+72406 60156 94 
6 145-34806 20139 04 
¥ —262-08507 46701 04 
8 423-48982 78720 85 
9 —612-78795 39166 28 
10 786:34127 71923 36 
11 —878-39229 00548 06 
12 832°55508 27743 51 
13 —649-04359 94610 06 
14 401-01836 81777 91 
ta —187-33116 59739 68 
16 61-77872 28094 79 
17 —12-76042 97736 49 
18 1-23748 46452 04 

Conclusion 


The polynomial coefficients given in the foregoing 
enable all nine of the functions 
It, th, fx, dx, Ko: >Ki, {x Ki dx, 
Ip—Ly, 1 -Li, fx (I) — Ly) dx 


to be solved to an accuracy of nine significant digits. 
The advantage over polynomial approximations 
published elsewhere [9] is the accuracy obtained; 


furthermore, the approximations for 
Jxtyds, fxKidx, h-Ly lly, [xth—Ly) dx 


are new. 

For practical applications, such as calculating the 
magnetic fields of transformers with concentric wind- 
ings, which have to be performed on a routine basis 
in the Company, as well as for axially symmetrical 
problems encountered in nuclear engineering [7], 
the accuracy afforded has proved a great relief to the 
calculators. 

(KME) A. KustTER 
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NON-LINEAR PLANE MAGNETIC FIELDS 


Numerous examples of linear magnetic fields are known. 
This article examines some of the non-linear, i.e. the saturated 
fields. For a certain class of magnetization curves some fields 
can be calculated exactly, as shown. 


Linear and Non-Linear Magnetic Fields 


AGNETIC FIELDS are said to be linear when 

the flux density B is proportional to the field 
strength H, that is to say when B = wH, the per- 
meability w being constant. This condition obtains 
when the magnetic field runs in air. If the lines of 
the field 
to be plane. Such magnetic fields have been tho- 


force run in parallel planes, is said 
roughly dealt with in the literature. A particularly 
suitable means of calculating them is the theory of 
functions. For experimental investigations the electro- 
lytic tank may be used. There is also an extensive 
literature on this subject. 

In contrast to fields in air, those in iron are subject 
to saturation. They are therefore non-linear and 
cannot be determined by the above methods. In the 
literature on the subject of saturation the main stress 
is laid on the magnetization curves. In this article, 
though, it is proposed to examine a few non-linear 
fields. 

These depend, not only on the geometrical form 
of their boundaries, but also on the magnetization 
curve. The present investigation is confined to a 
certain class of these curves, namely those conform- 
ing to the equation H = cB", where n is a fixed, 
real power. The measured curves approximate 
closely to this expression over a wide range. With 
this formula the problem is exactly soluble for a 
number of cases, such as axially symmetrical fields, 
for fields in a parallel strip, for the field in an angle, 
and for the plane dipole field. The calculated lines 


of force are illustrated graphically. 


538.122: 530,182 
The Magnetization Curves 


As mentioned above, the calculations are confined 
to the magnetization curve H = cB". By employing 
a suitable scale for H and B it is possible to make 
c= 1. This formula is chosen for the following 
reasons. If the flux density B is increased by a con- 
stant factor k, the field strength H, in the case of a 
single power, increases by the _ positionally in- 
dependent factor k". The lines of force and _poten- 
tial retain their shape. With other magnetization — 
curves, however, the formation of the lines of force — 
changes with the intensity. This limitation con- 
siderably simplifies the calculation in many cases. 
The field illustrated in the third chapter can be 
calculated for any desired curve, but this does not 
appear to apply to the other cases dealt with. 

The assumed expression also exhibits a remarkable 
symmetry. If we substitute 1/n for n, and H for B, the — 
equation of the magnetization curve is in effect 
unaltered. In the diagram the lines of flux and poten- 
tial are interchanged. This fact can be utilized to — 
obtain the potential lines from the formula giving — 
the flux lines. 

If the flux density B and the field strength H are 
both plotted to a logarithmic scale, the idealized 
magnetization curve becomes a straight line. Fig. | 
shows a measured curve. In the range from about | 1 to ~ 
20 kGs for the flux density, corresponding to an — 
interval from roughly 3 to 300 A/cm of the field — 
strength H, the curve is approximately straight. This 
is the most important range for practical applica- 
tions. Higher flux densities are seldom used and at 
low values the field strength is so weak that the 
absolute errors can be tolerated in the majority of — 
cases. This shows that the above formula is a suffi- — 
ciently good approximation for practical purposes. 
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Fig. 1. — Magnetization curve 


The flux density B and the field 
strength H are plotted to a logarith- 


[7 
ce 


mic scale. The curve is approxi- 


mately a straight line with the 
gradient 9. 


BROWN BOVER: 117228-"1 


Fig. 2. — Linear, axially symmetrical, spiral field, plotted 
‘for 2= 45° 


W — Flux lines 
@ = Potential lines 
r, p = Polar co-ordinates 
A= Angle between radius vector 
and field vector 


The powers n are obtained from the gradient of 
the straight lines. For the numerical calculation 
n = 9 was chosen. The adoption of n = | reverts to 


the linear case. 


Fig. 3. — Non-linear, axially symmetrical, spiral field 


Notation as for Fig. 2. BROWN BOVER! 117229+1 
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Axially Symmetrical Fields 


A field of this kind is composed of an eddy and 
a source. On account of the eddy, the magnetic 
potential increases by an amount H,, independent 
of the radius 7, in one revolution round the centre 
of rotation. Thus the azimuthal field strength is 
H,, = Ho/2zxr. If the source has an intensity of Bo, 
at the radius 7 the radial component of the flux 
density is given by B, = B,/2xr. 

If A is used to denote the angle between the radius 


vector and the vector of the field strength, H, 


Hsin A and B, = B cos d. Taking B and H from the 


no | 
= 
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above equations and substituting them in the for- 


mula for the magnetization curve, we obtain 


in which XK is a constant. In the linear case (n = 1), 
A is independent of r. The lines of force are logarith- 
mic spirals. In the non-linear magnetic field (n > 1), 
‘A, grows with r. For very small values of r the lines 
of force are almost radial. For larger values of r they 
are bent and tend to run almost azimuthally. The 
lines of force in the linear field are shown in Fig. 2, 


those of a non-linear field in Fig. 3. 


Fig. 4. — Linear field in a zone with parallel boundaries (parallel strip ) 


k = Parameter of the curve 
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Fig. 5. — Non-linear field in a zone with parallel boundaries 


The numbers included indicate the extraordinary increase in the potential ® with x. 


k = Parameter of the curve 


When general boundaries are imposed on the extent 
of the field, the magnetization curve and Maxwell’s 
equations yield non-linear partial differential equa- 
tions. In certain cases these may be converted to 


ordinary equations, two of which will now be discussed. 


The Field in a Zone with Parallel 
Boundaries (Parallel Strip) 


Let us consider a field which flows in one direction 


in one-half of a parallel strip, is bent back towards 


, 4 


the middle of the strip and flows back in the other 
half of the strip. In the linear case (Fig. 4) the equa- 
tion for a field curve may be expressed by x = 
k — log. cos y, where x and y are the co-ordinates 
along and across the strip, respectively; k is a con- 
stant. The width of the strip is made x. If one flux 
line is displaced parallel to the strip it coincides 
with another line. It can also be shown that the field 
exponentially with the co- 


strength increases 


ordinate x. 
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Fig. 6. — Linear field in an angle of 90° 


If we now postulate the hypothesis that the 
last-mentioned two properties also apply in the 
non-linear case, the partial differential equations 
are reduced to ordinary equations. It is apparent 
that no contradiction arises, and that the expression 
If we 


derivative dx/dy of a field curve, the following para- 


is thus justified. use p to denote the 


metric representation of a flux line is obtained 


= aresta Pa ze. 
y = arc tan p a 1 +R 
ip’ as on ae List p? 


The flux and potential lines are plotted in Fig. 5. It 
will be observed that the flux lines are compressed 
far less towards the edges of the strip than they are 
in the linear case. But in the middle of the strip 
they are bent much more. 


BROWN BOVERI 


117233+1 


Fig. 7. — Non-linear field in an angle of 90° 


The angular zone for a field running in saturated iron is marked 
by the hatched lines. 


The Field in an Angle 


The arrangement of the fields investigated in this 
case can be seen from Fig. 6 to 9. The angle | < « 
may be assumed to represent iron, the remainder air. : 
We are looking for a field which is completely in the — 
iron; hence the magnetic field in the air is ignored. 
Owing to the high permeability of the iron, except 
when heavily saturated (above about 18 kGs), this | 


does not cause any serious error. 


This case often occurs in practice in the vicinity © 
of edges, e.g. in transformers at the junction between — 
the core and the yoke and at the bends in the yoke, | 
in generators between the teeth and the yoke of the © 
stator laminations, or between the teeth and the 
rotor body. Particularly frequent are the right-angle i 
bends 2a = 90° and 2% = 270°. 

The formula obtained can be utilized to calculate 
the field of a plane dipole, by making « = —90° 
(Fig. 10 and 11). 
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Firstly the auxiliary quantity a is calculated from 


] 


Ries st ase 
(1 — 2a/7)2 


_ -(n+1) +) (n+ 1)? — 4nK 
ts 2K 


Fig. 11. — Linear dipole field 
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1. Along each ray passing through the centre 0 (i.e. 
along m = constant) the field strength increases 
as a power of the radius r (as, for instance, in 
Fig. 6). 


2. All flux lines are similar. 


It can be shown that when the magnetization curve 
H = B" chosen is used, both these properties are 
obtained; also in the non-linear case. With this 
hypothesis the partial differential equations are 
reduced to ordinary equations, which can be exactly 
integrated. 


The calculations can thus be reduced to the 
following final formulae: 


The flux lines can now be calculated in polar co- — 
ordinates r and m with the aid of the parameter 4 


from 
_ : a : lta . 
as A ag ™| 4 anal 
~ l atn ; lam ern 
r= aaa = Pav cot ) 2 (n+ a) 


The parameter A denotes the angle between the 


radius vector and the field line. The non-linear flux — 
and potential lines Y and 9, respectively, are illus- - 


trated in Fig. 7, 9 and 10, while Fig. 6, 8 and 11 
show the corresponding linear fields. Here too it — 


, 


may be seen that the flux lines are less severely com- 


pressed in the non-linear case. 

Although the selection of calculated non-linear 
field diagrams is still very modest, compared with | 
the linear cases, those reproduced in this article 
nevertheless indicate that it is not always impossible 
to calculate them. 


(KME) A. DALCHER 
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THE ANALOGUE COMPUTER, 
ITS EMPLOYMENT AND IMPORTANCE IN THE COMPUTER CENTRE 


The present article deals with the principle, operation, 
application and importance of the analogue computer in- 
stalled in the Brown Boveri Computer Centre. This machine 
is used for a wide variety of tasks, frequently acting as a 
simulator with normal equipment. With it the practical 
effect of any improvement made in the overall system 
can be observed directly. The article concludes with an 
example from the field of process control, pointing out the 
importance of the analogue computer as a means of predicting 


the operational behaviour. 


HE COMPANY has now possessed an analogue 
computer for about four years. This review of 
the experience gained in that period will give some 


idea of its application and importance. 


Principle of the Analogue Computer 


By way of introduction, the basic nature of the 
analogue computer may be briefly outlined, although 
there are plenty of references in the literature on this 
subject. It includes elements for all the elementary 
mathematical operations, such as algebraic sum- 
mation, multiplication, integration and the evalua- 
tion of non-linear functions. These operations can be 
represented with great accuracy in a frequency band 
between zero and about 50 c/s. 

The solution of a problem with the aid of the ana- 
logue computer may be divided into a number of 
phases, the first of which is an analytical phase—the 
installation or machine to be investigated being split 
into its elementary parts, although it must be added 
that such investigations are only undertaken for 
concrete objects. This analysis then provides a group 


681.14:530.17 


of mathematical equations or, in the particular case 
of a control problem, a block diagram containing 
all the important functions. Moreover, the analysis 
provides a closer insight into the nature of the system 
being investigated and, by omitting minor details, 
allows the principal functions to be recognized. 

The second phase is the synthesis in which, with 
the aid of the various elements of the analogue com- 
puter and according to the system of equations pre- 
viously determined, or the selected block diagram, 
the system is simulated, or a model analogue to the 
system under investigation is set up. From the ac- 
curacy with which the model is able to reproduce 
certain known functions of the object, it is possible 
to confirm the justification for the approximations 
chosen in the analytical phase. 

The next phase is the operation of the system, 
as simulated on the computer. The freely chosen 
parameters may be varied, readings are taken and 
any improvements in the general arrangement which 
appear necessary can be made. ‘The results so gained 
are then communicated to the interested departments. 

For all applications in the sphere of control engi- 
neering involving high accuracy, computing elements 
are required which possess still higher accuracy and 
stability than the elements being investigated. To 
cope with the steadily advancing perfection of 
engineering techniques, the analogue computer has 
also to be continuously extended, with the advan- 
tage that its capacity is increased at the same time. 
Practical experience shows that analogue computers 
can be employed to solve problems of increasing 
importance, necessitating the employment of a 


growing number of elements. 
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Importance of the Analogue Computer 
in an Industrial Undertaking 


The foregoing remarks dealt with the phases 
preceding the solution of a problem; these already 
provide quite a close insight into the inner structure 
of the task on hand. As soon as the circuitry and 
arrangement have been settled, the observation of 
the mutual effect of the various variables provides 
extremely valuable information, particularly as it 
is an easy matter to take readings simultaneously on 
several different elements. The effect of every 
improvement made to the arrangement can be 
checked immediately. 

In certain cases the circuitry on the computer can 
be set up to correspond exactly to the block diagram 
of the installation being investigated, without having 
The skilled 


operator can set up the analogue diagram with very 


to resort to mathematical methods. 


little extra effort, the great advantage of which 1s that 
the circuit is right before his eyes and he can correct 


any errors at once. 


Mathematical Theory and the 
Analogue Computer 


As far as the comparison between mathematical 
theory and analogue calculation is concerned, it 
may be established that the analogue computer is 
primarily an experimental device whose accuracy 
often has to be checked by mathematical methods. 
The computer may also be able to augment the 
theoretical approach to a solution, in that by vary- 
ing the variable parameters it provides indications of 


where the theoretical treatment has to commence. 


Analogue Computers as a Means of 
Instruction 


The ability to couple an analogue computer with 
normal equipment has already been mentioned. 
With this interesting combination it is possible to 
study the overall system under circumstances closely 
resembling those encountered in actual service. 
This leads to a further application of the ana- 


logue computer, namely as a means of instructing 


operating staff—a method which has been standard 
practice for some time in air transport undertakings, 
for example—especially when they have to control 
equipment which can be partly or wholly simulated 
on the computer. Particularly fruitful results have 
been obtained, for instance, in training staff con- 
cerned with the commissioning and operation of 
control systems, where it is necessary to demonstrate 
the manner in which an optimum control can be 
exerted, in the shortest possible time, without the 
installation suffering any ill-effects. For modern con- 
trol systems with an increasing number of stipulations, 
this method of instruction continues to gain in 
importance, especially as regards the time saved in 


searching for the optimum conditions. 


Operation of the Analogue Computer 


To judge the practical usefulness of the analogue 
computer the following aspects must be considered. 
As already mentioned, the computer is primarily an 
item of experimental equipment. Generally speaking, 
it is not suitable for routine work because the arrange- . 
ment of its elements is usually complicated and is 
subject to certain contingencies, even when operated 
by someone having perfect knowledge of the instal- 
lation. A certain time always has to be spent before 
it can be actually put into operation, owing to the 
number of checks which have to be performed. 
Since the structure of an analogue computer is very 
complex and sensitive, it requires a specialist to 
operate it. This operator then has the task of exe- 
cuting the more difficult problems and lending his 
assistance to other users of the computer. In this 
capacity he must be able to grasp a wide variety of 
technical problems in a very short time, and to 
analyse them, for which he must possess considerable 
mathematical ability. To enable him to deal with 
the often complex problems, he must also be able 
to think on synthetic lines. 

It should be most instructive for this specialist 
to be allowed to take part in the commissioning of 
equipment which he has previously investigated on 
the computer. It also enables him to widen his 
experience of the reciprocal effects of the various 


quantities involved. 
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Typical Application 


An example of the application of the computer is 
the investigation of a system for process control. For 
problems of this kind it is impossible to check the 
behaviour of the complete system in the laboratories 
of the firm supplying only separate elements of the 
system, because the controlled object, the main part of 
the whole system, is not available. With the analogue 
computer it is possible in most cases to simulate 
this object fairly closely, so that the dynamic response 
of the whole control system can be predicted with 
sufficient accuracy. 

In Fig. 1 the servomotor acts as the correcting 
unit and, via valve 2, varies the rate of flow. The 


servomotor is controlled by the dashpot regulator 
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(amplifier) 3. These parts are mounted on a test- 
bench in the workshop. The other parts of the control 
system, i.e. the electronic amplifier 4 and the ana- 
logue computer 5 are situated in the Computer 
Centre about 100 m away from the test bench. As 
indicated in the sketch, a multicore cable acts as 
the connecting link between the two parts of the 
control system. 

In the computer the amplifier A, and the connec- 
ted resistors and capacitance simulate the controlled 
object, which in this case can be represented by a 
single time constant. A, and A, represent the transfer 
functions of the measuring system, which can be 
described by two time constants. A; produces the 


reference input depicted in Fig. 2, while A, acts as 


summation element. A, merely serves a purpose 
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Fig. 1. — The control circuit being investigated 


1 = Servomotor 

2= Valve 

3 = Dashpot regulator (amplifier) 
4 = Electronic controller 

5 = Analogue computer 

A,-A, = Amplifiers 


valve of 1 


M,—M, = Setting potentiometers 
I,-I, = Control currents of 3 
I, communicates the feedback of the control 


I, communicates the feedback of valve 2 
B,—B, = Control solenoids 


The components | and 3 are situated 

in a part of the factory, while 4 and 

5 are in the Computer Centre. The 

signals are transmitted over a multi- 
core cable. 
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Fig. 2. — Traces obtained during an extensive series of tests 


a: Desired value 
b: Actual value 
c: Difference between a and b 
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associated with the instrumentation and produces 
the difference between the reference value and the 
actual value. 

The upper part of Fig. 2 shows the master or 
desired value, the middle section the actual value, 
and the lower section the difference between the two. 
It will be observed that the actual value responds 
to changes in the desired value with a considerable 
time-lag, but that there is very little overswing. 

As a result of this investigation it was possible to 
predict that the components of the control system to 


be delivered would exhibit the correct response. 


(KME) P. CoROLLER 
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METHODS OF INVESTIGATING ELECTRICITY NETWORKS 


After briefly mentioning two mathematical methods, the 
author describes the investigation of electricity systems with 
the a.c. network analyser, the remarks being illustrated by 
reference to practical examples. Some useful tips are given 
regarding simulation techniques. These are followed by 
reference to some more recent methods, e.g. coupling the 
network analyser to the analogue computer, and digital 
computation of the network by means of the program- 
controlled computer. 


On the Principles of Investigating 
Electricity Networks 


NYONE confronted with the task of calculating 
the conditions in an electricity network will 
immediately recognize the characteristic difficulty in 
the handling of such problems: Having regard to the 
particular configuration, every network requires a 
fresh calculation, starting from scratch, involving an 
extremely large number of variable quantities. 

For the easier investigations two possible methods 

immediately spring to mind. The first method is to 
simplify the configuration of the network by elimi- 
nating the junction points. The original network is 
replaced by an equivalent network, enabling the 
problem to be solved more easily. It is also possible to 
carry out a reverse transformation, thus determining 
the voltages and currents in the original network 
step by step [1]. 
_ The second method is to set up the linear system of 
equations from which the currents or voltages of the 
various links forming the network can be determined 
[2]. 

However, for the following reasons it is impossible 
to employ such a method when calculating by hand 
or with an ordinary office calculating machine, except 
perhaps in certain exceptional cases. 

The given quantities are usually not the resistances, 
voltages or currents, but mixtures of resistance values, 
absolute values of voltages and some of the powers. 
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The unknowns are then the remaining voltages and 
powers. 
Furthermore, as regards the electricity networks, 


the following development may be observed: 


| 


The networks are continually growing; 
— Interconnection is increasing; 


— The number of points where junctions are made 


with neighbouring networks is growing; 


— The number of infeed points is increasing, either 
from the network’s own generating stations or from 


neighbouring networks; 


— The number of consumers or shunt links in the 
interconnected system of networks is always being 
enlarged. With the use of progressively higher 
voltages, the importance of capacitive shunt bran- 
ches grows too, with the result that the number 


of links in the system continues to increase. 


Hence, as a rule, it is only possible to investigate 
electricity systems with the aid of analysers or digital 


computers. 


Analogue Devices 


The d.c. network analyser, the analogue computer 
[3] and the a.c. network analyser [4] are all analogue 
devices. The first two operate in principle with direct 
current, the third with alternating current. W. J. 
Karplus [5] places the network analysers in the 
category of direct analogue devices, but the analogue 
computer in the indirect category. This classification 
expresses the main difference in their operating 
principle. 

Direct Analogy 
In a direct analogy one and only one physical 


quantity in the simulation corresponds to one phy- 
sical quantity in the original. Thus, for instance, in 


the thermo-electric analogy, temperature difference 
corresponds to voltage difference, thermal resistance 
to ohmic resistance and the flow of heat to current. 
When determining the natural frequency of a shaft 
subjected to pulsating torques, the moment of inertia 
corresponds to the inductance, the angle to the 
current, and the reciprocal elastic torque to the 
capacitance. If a dual image of an electricity system 
is set up, a current source corresponds to a voltage 
source, an inductance to a capacitance and vice 
versa. The particularly lucid case, in which corres- 
ponding values are represented by the same physical 
quantities in both original and image, are contained 


as a special case in this direct analogy. 


Indirect Analogy 


In contrast, in the indirect analogy, different 
quantities in the original are represented by a single 
physical quantity in the set-up. Thus, for instance, 
when setting up an electricity system on the analogue 
computer, the current and its integral, the quantity 


of electricity, are determined by measuring a voltage. 


The next chapter will describe the investigation of 


electricity networks, using an a.c. network analyser. 


The A.C. Network Analyser 


The Brown Boveri a.c. network analyser (see Fig. 3 
on page 285) contains several voltage sources, auto- 
transformers and variable impedances, which are 
used to simulate the generators, transformers, trans- 
mission lines, shunt reactors, synchronous capacitors, 
capacitances and loads in the network being investi- 
gated. It is equipped with metered and unmetered 
jumpers, as well as shorting links. The available ele- 
ments are listed in the table of active and passive 


elements on page 315 opposite. 


Of the twelve generators, five exhibit a reactance 
in the range from 0 to 6-21 per unit, ten generator 
elements can be loaded up to 10 per unit and two 
up to 20 per unit. Six capacitive elements can be set 
at up to 7-21 per unit. Other elements which can be 
plugged in at the rear of the analyser allow the 
maximum value to be increased to 11-21 p.u. A num- 


ber of the auto-transformers allocated to the load 
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circuits can be varied down to 0-4 p.u. Furthermore, 
these auto-transformers and those in the succeeding 
lines can have their voltages varied in steps of 0-005 
p-u. The number of insulating transformers has 


increased to ten. 


Having made all the necessary settings, and on 
completion of balancing, minor changes in the cir- 
cuitry (e.g. opening lines, switching on loads, caus- 
ing short circuits, etc.) and any adjustments to the 
amplitude and phase angle of the voltage sources 
can be made from the control desk. Special measuring 
amplifiers in the generator elements and in the 
metering desk prevent the instruments from affec- 
ting the image of the network. All passive elements 
can be adjusted in decades. As a direct analogy of the 
original network, the analyser is an extremely adap- 
table device and provides a very clear picture of the 


network being investigated. 


Reference to Completed Investigations 


As examples of the numerous investigations which 
have been carried out with the a.c. network analyser, 


the following are worth special mention. 


Flow of power 


Interconnection of the networks of a number of 
electricity undertakings under normal operating 
conditions and with special connections. The investi- 
gation of successive stages in the expansion of net- 
works at home and abroad, including railway trac- 
tion networks, urban supply systems which may al- 
ready be in existence, or due for reorganization or 


extension, and finally factory supply systems. 


Faults 


Included under this heading are investigations to 
determine the short-circuit powers in the event of 
three-phase short circuits (e.g. to determine the cir- 
cuit-breaker rating needed). The a.c. network ana- 
lyser can be used to represent networks with negative 
reactances (three-winding transformers) and bran- 
ches of any desired impedance angle. The problem 
which had to be investigated most often was the 
single-pole earth fault in a network with earthed neu- 
tral. Other investigations were concerned with the 
line-to-line fault, single-phase open circuits, asym- 


1 
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Active and passive elements of the a.c. network analyser 


Total R St 
No Designation uan- on Messen ore a 
} q : which be set Remarks 
tity in the per unit system 
1 | Generator element 12 Voltage U 0-2°5 continuous | Output terminals shorted 
Angle p 0-360° when “Off” 
2 | Generator reactance} 12 7 Reactance wL 0-1-220 0-001 — 
0-001 from 0-1-220 
5 0-6:21 
0-01 from 1-00-6-21 
3 | Line element 114 90 Resistance R and 0-1-220 0-001 _ 
0-001 from 0-1-220 
4 24 reactance wL 0-6-21 
0-01 from 1-00-6-21 
5 | Capacitive 60 12 Susceptance Y 0-0-621 0-001 
element 42 0-1-221 ca 
0-002 from  0-1-220 
6 0-7-21 
0-01 from 1-00—7-21 
7 | Load element 27 Resistance R and 0-2-16:2 0-01 ~ | Connection in series or 
reactance wl parallel 
8 | Auto-transformer 12 0-7-1°31 0-01 
for load element 27 Transformation ratio with 2 “Off? positions 
15 0-4-1-315 0-005 
9 | Auto-transformer 12 Transformation ratio 0-7-1-315 0-005 = 
10 | Insulating transfor- 10 6 Ratio 1:1 
mer 2 ae — — Leva 
2 1:V3 
11 | Unmetered jumper 48 — — oo non-measurable 
12 | Metered jumper 24 — — —_ measurable 


For use in the analyser the per unit system is based on the following rules: The rated value (15 V) of the output voltage of the 
generator element is denoted as 1, the rated value of the current (30 mA) also as 1. All passive elements contained in the analyser 
are appropriately calibrated according to the per unit system and marked. The unit of resistance is 15/003 = 500 ohms and that of 
the conductance 0-03/15 = 0-02 mho. 


J 
_ metric branches connected in series (a capacitor 


short-circuited in one phase) and faults in star/delta 
transformers. These investigations are carried out 
th the methods employing symmetrical components 


[6] or «, 6B, 0 components [7]. When the symmetrical 


components are employed, the generator reactances 
in the positive, negative and zero-sequence systems 


can be simulated very exactly. As a normal rule, the 
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negative and zero-sequence systems only contain 
passive elements and can therefore be replaced by 
pre-determined branches or polygonal networks, 
depending on the case in point. In contrast, active 
elements have to be inserted in both the « and the f 
system. However, the «, 8, 0 components, compared 
with the symmetrical components, offer the advan- 
tage that, for the equivalent circuits for any kind of 
multiple fault, only insulating transformers are requi- 
red and no means of rotating the voltage and current 
vectors [8]. A further problem belonging to this 
category was that of resonance with harmonics of the 
5th, 7th, 11th and 13th order. Long lines in the 
original network were represented by 7 lines. The 
network was simulated for the 6thand 12th harmonics 
and the frequency varied continuously between the 
5th and 7th, or between the 11th and 13th order. 


Stability 


Particularly interesting problems are encountered 
in stability investigations. ‘The step-by-step method 
demands great concentration, but it underlines the 


special advantages of the analogue method. 


By continuously determining and adding up the 
angular differences the engineer is simultaneously 
able to realize why individual machines and the 
whole system behave the way they do and not in 
some other manner. He gains insight into the effect 
of the various parameters and, by different measures, 
is able to control the process so that the desired effect 


is achieved. 


The problems dealt with have concerned the sta- 
bility of networks with synchronous generators in 
hydro-electric, steam and diesel-powered generating 
stations, with synchronous motors, synchronous capa- 
citors and frequency converters; also with symmetric 
and asymmetric faults with progressive adjustment of 
the loads, as well as cases which could only be simulat- 


ed artificially (dummy loads on the generators). 


Special problems 


These include numerous transfigurations of net- 
works into branches or equivalent triangles or poly- 
gons; also compensation circuits, dimensioning the 
means of modulating for variable frequency, measur- 


ing equipment for unbalanced loading, circuits for 


converting the component values into phase values 
(symmetric components) and natural frequencies of 


shafts subjected to pulsating torques. 


The smallest network investigated by the analyser 
contained 5 branches, the largest over 200. The statis- 
tical break-down of the network investigations is as 


follows: 


57% for the determination of power flow, and the 


distribution of voltage and current, 


25° for the investigation of symmetric or asymme- 


tric faults, 


10% for investigations of steady-state and transient 


stability, 


8°4, for electrical and mechanical problems of a 


special nature. 


Special Artifices 


Three-winding transformers 


These can be simulated by a reactance star [9], 
in which it is by no means rare for negative reactances 
to occur. A small negative reactance corresponds to 
a large capacitance, which sometimes exceeds the 


maximum value available on the analyser. 


By utilizing a star/delta transformation a reactance 
star can be converted into an equivalent delta. If 
Ys1, Ys2, Ys3 denote the susceptances of the limbs of 
the star, and if the subscript 3 is allotted to the limb 
of the triangle across limbs | and 2 of the star, its 
susceptance Yp3 is given by 


Ys, - Yso 
Ys, + Yso+Ys3 


Yp3 = 


Multiplying both sides of the above equation by Ys3 
yields 
pwr Vie Vos say 
Vide Gi ee lee nae 
Vsti sales 


Thus if the susceptance Ys3 of one limb of the 
star is large, from the absolute aspect, that of the 
opposite side of the triangle Yp3 will be correspond- 
ingly small, i.e. more favourable for the analyser. In — 
that case it is preferable to replace the reactance star 
with a negative reactance by the equivalent delta. 


; 
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Transformers with combined direct and quadrature regulation 


In this form of regulation the direct-axis com- 
ponents are obtained by proportionately raising or 
lowering the transformation ratio by means of an 
auto-transformer. The quadrature-axis components 
can be realized by simulating the whole system twice, 
once in an « system and once in a f system [10]. With 
an angular displacement of, say 60°, the ratio of the 
direct to the quadrature-axis component is b:/3. 
The components are therefore conveniently com- 


posed with insulating transformers having the ratio 


1:/3. 


Leg of a polygon with negative resistance! 


The impedances of individual legs of a polygon 
equivalent to a network sometimes contain real com- 
ponents with a negative sign. This state of affairs 
occasionally occurs during the transfiguration into 
an equivalent delta of several consecutive line sec- 
tions with infeed reactances. If no special elements 
are available for simulating the negative resistances, 
it is possible to circumvent the negative resistance of 


the leg of the polygon by using passive elements. 


In the triangle (Fig. 1) R,, X, and R are all positive 
values. The leg AB has an impedance with a negative 
real part, while AC and BC have the same ohmic 
resistance. This delta can be converted into a star, 
whose impedance only contains positive real parts, 
provided the inequality 


G)<*<3 +(e) 


is fulfilled. 


(1) 


Any side of an n-sided polygon forms triangles 
with the (n—2) points not on the side in question. If 
one of these sides has an impedance with a negative 
real part, of the form Z = —R, + jX;, first a resis- 
tance R is chosen, which fulfils the above inequality 
(1). Then from the (n—2) triangles a suitable tri- 
angle is sought, whose legs connecting to an opposite 
point of the polygon fulfil the following condition: 
The real parts of the conductances of the two legs must 
not be less than the reciprocal of R. The resistances 


1 We are indebted to W. F. Grauam, Bureau of Reclamations, 
Denver, USA, for bringing this method to our attention. 


R can now be split off these two legs in the form 
of parallel resistances. These, together with the leg 
having the impedance Z, form a triangle which can 
be transformed into a star having no negative resis- 
tances, so that now all legs of the new equivalent part 
of the network can be simulated by passive elements. 

By way of example, let us assume a triangle whose 


legs have the following impedances, as shown in 
Fig. 2a: 


Z,= —R, +jX,=-01 + j3-6 
Z,= 0-06 + j0-32 
Z,= 0-08 + j2-1 


To determine the resistances R belonging to the 
first leg AB, the inequality (1) must first be set up. 
With (R,/2) = 0-05 and (X,/R,) = 36 it reads 


0-05 < R < 64-85 
and correspondingly, for the conductance G 


0-01542 < G < 20 (2) 


The susceptances of the legs AC and BC work out 
to 0-56604 — j3-01887 and 0-018114 — j0-47550, re- 
spectively. Their two real parts are larger than the 
lower limit in (2) so that the negative resistance of 
the first leg can be eliminated. If G is made 0-018114, 
the legs AC’ and BC can be resolved into two com- 
ponents each, namely Z;, R, and Z/, R;, as shown 
in Fig. 2b. The dotted triangle consisting of the leg 
with the negative impedance Z, and the two ohmic 
resistances R, = R, is now transformed into the 
equivalent star. In this way Fig. 2a is converted into 
the network in Fig. 2c containing no negative 


resistance. 


A B 
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Fig. 1. — Delta network 


The leg AB contains a negative resistance. The other legs 


have equal ohmic resistances. 
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Fig. 2. — Determination of an equivalent network containing no 
negative resistances 
Values expressed in the per unit system 
a: Original network. Leg AB contains a negative resistance 
b: Resolving the legs AC and BC 


c: The resultant equivalent network obtained by the delta-star 
transformation of the dotted triangle 


A.C. Network Analyser in Conjunction 
with the Analogue Computer 


The analogue computer is now accommodated in 
the room next to the network analyser. Preparations 
are at present being made for coupling these two 
machines together. When the analyser operates alone, 
the amplifiers representing generator elements are 
fed from the RC oscillator, the amplitude and phase 
angle of the voltage being adjusted at the source or 
from the desk. When the two machines operate 
jointly, these settings on the desk are blocked. The 


voltages required are now fed to the individual ampli- 
fiers of the analogue computer via special components, 
their amplitude and phase angle being set on the 
analogue computer. 

The combination of the network analyser and the 
analogue computer opens up interesting possibilities 
for the investigation of electricity networks or other 
systems. The computer now has a network available, 
so to speak. The analyser, which was mainly used for 
steady-state tasks, can now be augmented by auto- 
matically controlled machines. Thus, in the near 


future the following problems can be solved: 
— Problems associated with automatic control, in- 
cluding the repercussions from the network; 


— The combined effect of control processes consider- 
ing several machines connected to the same net- 


work; 


— The effect of switching large machines on to a 


network ; 


— The effect of the control system on the magnitude 


of the short-circuit currents; 


— The curve of the short-circuit currents with respect 


to time; 


— Investigation of the transient stability of a number 
of machines, the oscillation curves being plotted _ 


automatically ; 


— The effect of voltage control on the transient 


stability ; 


— Repercussions of its control system on the prime 


mover. 


Digital Methods 


The solution of the problems mentioned in the 


opening chapter already involved digital methods. _ 
Although quite neat methods have been devised for 
network problems solved by means of desk-type cal- 
culating machines [2], their use is usually confined 
to very small networks. When a program-controlled 
digital computer is employed, digital methods are 
also suitable for large networks. Thus a new, very 


versatile instrument becomes available for the investi- 


gation of electricity networks. 
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This computer is ideal for calculations which recur 
cyclically with a certain regularity. In this way, for 
instance, shunt branches of consecutive line sections 
can be eliminated by the star-delta conversion of a 
group at a time, the network being then systemati- 
cally investigated with the aid of matrices [11]. The 
following method, of which only a brief description 
is given, is useful for investigating the flow of power 
and stability; this is the iteration method utilizing 
the equations of the junctions [12, 13]. 


Power Flow 


Starting with the estimation of the voltages, the 
currents flowing into each leg of the network from a 
junction point are calculated and added together. 
The state at this point may be either of two possi- 


bilities: 


1. The junction is unloaded, so that the sum of the 
line and capacitance currents must be zero. If this 
condition is not fulfilled, the voltage has to be 
changed. The condition of zero current provides 
two equations from which the supplementary 
voltage can be determined. 


2. The particular junction is connected to a gener- 
ator or load. In this case there are two more 


alternatives: 


— The incoming (or outgoing) active and/or re- 


active power is known, or 


— The active power and the absolute value of the 


voltage are given. 


Here, as with the unloaded junction, two equations 
can be established, enabling the conditions to be 
fulfilled when the voltage at the junction has been 
corrected. 

Each junction in the network is then dealt with 
in this manner, the voltage being successively cor- 
rected and this correction stored if it exceeds the 
correction of the preceding point. Having dealt with 
all junctions, the next question to be settled is whether 
the maximum correction is smaller than a certain 
limit. If this is so, a stop command is issued; if not, 
a new cycle commences, the procedure being con- 
tinued until the correction is sufficiently small. Fol- 
lowing the calculation of the voltage, the flow of 
power is determined. The result appears in the form 


of a printed table, the amplitude and phase angle 
of the voltage appearing as a per unit value for each 
junction, while the powers of the feeding generators 
and the loads, as well as the powers flowing towards 
other junctions, are expressed in MW and MVar. 
In addition, the power losses in the network can also 


be printed out. 


Stability 


To determine the oscillation curve of synchronous 
machines it is necessary to know the amplitude and 
phase angle of their e.m-f. in the steady state, i.e. 
before the occurrence of a fault. Following the deter- 
mination of the generator outputs, by the analogue 
or digital method, these e.m.f. are easily determined, 
provided the transient reactances are known. The 
solution is then obtained by integration of the differ- 
ential equations of the synchronous machines, the 
active powers of the generators being determined 
each time by transferring to a power-flow program. 
This procedure may at first appear rather clumsy. 
At the beginning of the initial time phases several 
cycles are necessary owing to the changes in voltage. 
But subsequently the voltage differences are quite 
small in view of the short time interval. The junction 
voltages determined at the beginning of an interval 
provide a good approximation for the subsequent 
time value, so that, usually, only a small number of 
cycles are required. The line affected by the fault 
can be disconnected or reconnected, as desired. 

This method, however, offers the advantage that 
the network is available in its original state and not 
in a simulated state. This makes it possible for any 
supplementary values which may perhaps be re- 
quired to be determined, such as the voltages at 
individual junctions or the currents in circuit- 
breakers. 

The tape on which the results are printed out gives 
the angles of the e.m.f. and their derivatives relative 
to any chosen machine, or with reference to a fixed 


voltage (infinite supply assumed). 


Analogue and Digital Methods 


Both methods—the analogue and the digital— 
exhibit some striking advantages as regards their 


application for the investigation of networks. 
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With the analogue method every parameter, i.e. 
every load, every line, every transformer and every 
generating station is available for carrying out a 
particular adjustment. If the results of a simulation 
are presented on paper, questions such as the follow- 


ing may arise: 

What happens when 
— this transformer capacity is doubled ? 
— this line is only run as a single circuit? 


— these two sets of busbar systems are operated 


separately ? 


By making the appropriate settings on the ele- 
ments (3-6 selector switches have to be operated on 
the analyser) or by forming two junctions, a result 
is obtained which shows whether this change serves 
a useful purpose, or whether its effect is so small as 
to be insignificant. Especially during investigations 
connected with the extension of networks in stages, 
the analyser is an ideal aid to the supervising engineer 
as it allows him to look for the most favourable solu- 
tion. He knows exactly what means he can call on in 
the particular year and he can utilize these in the 
most suitable combination. All settings on the active 
and passive elements, as well as the circuit, can be 
quickly and reliably checked. The measured power 
values are entered in on the network diagram direct 
and immediately checked by adding round the junc- 
tion. This representation gives a very clear picture of 
the flow of power in the network. 

The digital method allows large, completely 
meshed networks with numerous feed points and 
loads to be investigated. A further advantage is the 
higher accuracy attainable with this method, al- 
though this causes a corresponding increase in the 
time taken (more cycles are needed). It is thus easier 
to compare different variants with one another as 
regards their respective losses. 

A feature of the digital method is that it allows the 
methods used previously to be refined by omitting 
simplifying assumptions, and the process to be spread 
over a longer period. Furthermore it will allow prob- 
lems to be investigated which presented insurmount- 
able difficulties previously, either owing to the time 
taken or the size of the network in question. Com- 
pared with the analogue method, the digital method 


appears to be most suitable for those networks for 
which a preset plan can be given, according to which 


the investigation has to be carried out. Having re- — 


gard to the wide variety of tasks, which the arrange- 
ment, operation and extension of electricity networks 
impose, it is advantageous to be able to adopt the 
one approach or the other, with the aid of the ana- 
logue or the digital methods. 


(KME) E. FAssLer 
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THE USE OF DIGITAL COMPUTERS 
TO EVALUATE MEASUREMENTS ON THERMAL TURBO-MACHINES 


To evaluate measurements, the same calculation has to 
be repeated for each measurement, merely using other 
numerical values. If this work is performed by an electronic 
computer, the main advantage is that the results are available 
within about one hour of the readings being taken. Conse- 
quently better use can be made of the test equipment. 
Depending on the nature of the machine being investi- 
gated, the computer program can be made up of a number 
of sub-routines. The formation of the various sub-routines is 
described in the present article. 


HEN THERMAL turbo-machines are being 
investigated, measurements are taken of quan- 
tities, from which it is not generally possible to ob- 
tain a direct assessment of the machine. This involves 
a number of calculations, normally referred to as 


“evaluation”’. 


Object of the Calculations 


These calculations should provide the results for 
the following problems: 


1. Many of the interesting quantities cannot be read 


direct, but must be determined from other 
measurements, i.e. by an indirect method. Thus, 
for instance, to measure the mass flow of liquids 
or gases the pressure difference at a throttling 
device, and the pressure and temperature before 
the device are measured. The mass flow can 
then be calculated from these readings and co- 


efficients determined by calibration tests. 


2. Other quantities, by virtue of their definition, 
can only be determined mathematically from 
other measurements. For example, the energy 
input to the gas in a compressor, or the energy 
extracted in a turbine have to be determined. 
This necessitates the observation of measurements 
at a number of points—at least two—on the ma- 


681.14: 621-135 


chine, from which the change in the state of the 
gas can be followed. Then, by allowing for certain 
properties of the gas, the difference in energy can 


be calculated from the change in state. 


3. From the individual quantities certain ratios or 
dimensionless quantities have to be derived. Thus, 
for instance, the efficiency is the ratio of the output 
power to the input power; or the load on the 
machine is measured by relating the energy con- 
verted in it to an energy referred to the peri- 
pheral speed of the rotor. The results of tests can 
be expressed more simply in this manner and are 
more generally valid, because they apply directly 
to “‘similar’’ machines. In this way the assessment 


of the investigated machine is simplified. 


Advantages of Calculation by Machine 


The calculations necessary for the evaluation of 
measurements, depending on whether the dimension- 
less quantities for the whole machine are required 
or whether they are also required for individual 
parts of the machine, can be very extensive. In con- 
formity with the number of runs which have to be 
performed by the particular machine, these calcula- 
tions have to be constantly repeated, using the same 
procedure but different numerical values. By em- 
ploying a digital computer the following advantages 
are gained over calculating methods of a non- 


mechanized nature: 


— The results of experiments are available much 
more quickly than when the measurements are 
evaluated ‘“‘by hand’. This is very important, 
because the test machine cannot be modified for 
a further series of measurements until the results 


of the first series have been evaluated. 
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— The accuracy of calculation is increased. 
— The reliability of the calculation is improved. 


— The cost of calculation is reduced, provided the cost 
of programming can be spread over a sufficiently 


large number of calculations. 


— The staff liberated by changing over to mechanized 
computation can be employed for other, more 


interesting duties. 


These considerations resulted in a change being 
made to the evaluation of experimental results by 
computer as long ago as the spring of 1957. To begin 
with the calculations were performed by an IBM 650 
computer of the International Business Machines 
Centre, Zurich, later by the Siemens digital com- 
puter type 2002 installed in the Brown Boveri works 
in Baden. 


Typical Application 


A typical calculation of this kind is the evaluation 
On_ these 


machines the effect of a large number of geometrical 


of measurements on model machines. 


quantities on the characteristics of the machine, 1.e. 
turbine or compressor, have to be determined. In 
view of the large number of experiments which 
have to be performed, the advantages of mechanizing 
calculation become particularly pronounced. 

The principle is as follows: During the experi- 
ments, the readings are recorded on a report sheet. 
Immediately the experiments are complete, an input 
punched tape containing all the readings has to be 
made, which together with the tape giving the instruc- 
tions for computation is read into the computer. 
On completion of the computation, a punched tape 
with the results is turned out, the actual table of 
results appearing in clear on the sheet printer. The 
time taken by the machine for calculation, including 
the input and output times, is about 10 minutes for 
15 measuring points. A skilled engineer, calculating 
by hand, would require about two days for this work. 
Thus the interval between successive series of experi- 
ments can be reduced by two days. As a result of 
using the computer, the test equipment is better 
utilized, because far more experiments can now be 


performed in the same space of time. Moreover, 


according to a rough estimate, the cost of calculating — 


by the computer is only about half the cost of calcu- 
lating “by hand’’. 

It is hoped that, at some later date, it will be 
possible to eliminate one of the intermediate steps 
described, when all measurements, instead of making 
the detour via the report sheet, are punched direct 


on the tape by a punch controlled by the measuring 


elements. Parallel with punching, the results would 


then have to be printed out, thereby providing a 


visual check. 


Composition of the Computer Program 
from Sub-Routines 


How must the computer program be organized — 


in order that the above advantages may be enjoyed ? 


Above all, it must be borne in mind that the experi- — 


ments have to be performed on machines of very — 


different kinds, i.e. turbines and compressors, which 


may be of axial or centrifugal design. For this reason ~ 


alone, the arrangement of the measuring points and 


the procedure for evaluation is bound to be different. 


Furthermore, many of the experiments differ in their — 


prime purpose: For instance, for guarantee tests 
only the main characteristic data are determined, 
such as the overall efficiency, and the throughput 
as a function of the load, while for tests in the develop- 
ment stage greater attention has to be paid to indi- 
vidual measurements. 

One method is to prepare a separate program for 
each case, intended only for a definite machine and 
a definite purpose. In this case the cost of preparing 
the program and checking it are relatively heavy, 
especially when only a small number of measure- 
ments are involved. On the other hand the program 


can be organized in such a manner that it takes all 


possible conditions into account. But such programs | 
become very large and unwieldy. For this reason — 
it has proved preferable to collect identical parts of — 
the program in sub-routines, which are already pro-_ 


grammed and checked. They are stored in the library 


in the form of packs of punched cards, or as punched — 
tapes, and can be combined to form a program adap- | 
ted to the particular test object and experimental 


set-up. 


MAY/JUNE 1961 THE Brown Boveri REVIEW 323 


Available Sub-Routines 


The parts of the computer program suitable for 
combination as sub-routines will now be described. 


1. Calculation of the mass flow of water or air from 
measurements obtained with the standard nozzle 
in accordance with the procedure laid down by 
VDI [1]. The formula given in these rules 


M = aAgke/20, Ap 


applies to incompressible fluids and compressible 
gases. In it the following notation is used! (see also 
Fig. 1): 


M=mass flow 
«= coefficient of discharge* (including correction 
for roughness of the pipe-wall) 
Ap =area of the nozzle at ambient temperature 
k =correction for temperature 
€ = expansion factor 
0, = density of the fluid before the nozzle 


Ap= pressure difference at the nozzle 


The formula can be broken down into a factor K, 
which must be calculated alike for incompressible 
liquids and compressible gases 


K = K (d,D,4,t,) 


where d= diameter of the nozzle 
D = diameter of the pipe 
t, = temperature before the nozzle 
2 = coefficient of superficial expansion 
(= twice the coefficient of linear 


expansion) 


and a factor B, which expresses the nature of the 
fluid whose quantity is to be measured. 


Incompressible B;, = B;, (Ap, t,) 
Compressible B, =B, (Ap, t,, p,, x, X, D, d) 


in which x is the ratio of the specific heats at con- 
stant pressure and at constant volume 
X is the absolute humidity. 


Then M=KB 


1 All the notation used in the article is summarized at the 
end of the article, see page 329. 

* This coefficient, as defined by VDI, is not identical with 
that defined according to American or British Standards, because 
it includes the velocity of approach factor. 


Ap 
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pi 


Fig. 1. — Arrangement of standard nozzle in the test pipe with 
the measuring points for determining the mass flow 


For notation, see text. 


Since the pressure before the nozzle is only used to 
calculate the quantity of compressible gases, it can 
be arbitrarily made equal to zero for incompressible 
liquids, and utilized for controlling the computer 
program (Fig. 2). 
The functions determined experimentally 
k =k (é, A) 


Ow ee ; 
«=«[(p), >| 
a id 
n) (>) 


where p, = p, — Ap, the pressure after the nozzle 


& 

I 

ie») 
<i 
oe a 
>| 
La | 


and oy = the density of the water 


Be\Ap, t), p,, x , X, d, D) | 
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Fig. 2. — Flow diagram for computation of mass flow, sub- 
routine I 


For explanation, see text, 
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are represented in the program by approximation 
functions. This point will be dealt with in greater 
detail in one of the later chapters. 

The flow diagram (Fig. 2) shows the computation 
procedure with the inputs and outputs. For the 
input data, those shown above the dotted line are 
characteristic of the measuring set-up and remain 
constant for a particular arrangement. Below the 
dotted line are the measurements which are different 


for every measuring point. 


2. To calculate the state of air or combustion gases 
from a measured temperature and a measured 
pressure. It has already been mentioned that to 
calculate the energy put into or taken from a gas, 
the state of the gas must be known in at least two 
positions on the machine. By measuring two pro- 
perties of state, e.g. the pressure and the temperature, 
the state of the gas at the particular point is un- 
ambiguously determined. But often the gas velocity 
at the point is so high that detecting elements 
mounted at a definite point record different values 
(owing to the flow in front of the measuring element 
being retarded and compressed) to those which 
would be recorded by detectors entrained by the 
flow. Thus the true values have to be obtained by 
correcting the measured values, using correction 
factors which indicate how much of the dynamic 
head is included in the measurement. The correction 
factor is unity when the gas is retarded to stagnation. 

For the temperature sensing element this is known 


as the dynamic correction factor 


3 he ala 


o= eT 


where 7J'= temperature of the gas 
T,, = measured temperature 


7* = temperature at the stagnation state 


and, accordingly, for pressure measurements, the 
dynamic correction factor of the pressure sensing 


element 


o= Pm —P 

Me 
with corresponding meanings for the various sym- 
bols. These correction factors have to be determined 
by calibration tests and may be regarded as constant 


within a certain speed range. 


To correct the measurements it is also necessary 
to know the velocity at the measuring point and the 
enthalpy difference Ah corresponding to the kinetic 
energy, in order to be able to calculate the differences 
T* — T and p* — p.2 To accurately determine the 
velocity at the measuring point, at least one further 
quantity, e.g. the stagnation pressure, would have to 
be measured. Often it is quite impossible to install all 
the necessary detecting elements at one point or at 
comparable points, and in any case they complicate 
the layout of measuring equipment. Consequently, 
for the sake of the correction of the measured 
quantities, the velocity is often calculated approxi- 
mately from the measured mass flow M, the cross- 
sectional area A, and the state of the gas, provided 
the angle of flow « can be taken as known. 

Hence this calculation requires those properties of 
state which have to be computed with the aid of the 
program. For the calculation by the computer it is 
most convenient to adopt an iteration method 
(Fig. 3). The measured quantities are regarded as a 
first approximation to the true property of state. 
With these an approximate value is calculated for 
the velocity at the measuring point and a first cor- 
rection made to the measured values, which then 
act as the starting point for the next approximation. 
The iteration is discontinued when the improvement 
in the enthalpy difference Ah is less than a definite 
value ¢ from one step to the next. 

Again in the flow diagram the input values given 
by the geometry of the measuring point and by the 
sensing element are shown separate from the measure- 
ments. With the latter the gas constant R must also 
be included, the value of which varies with the 
composition of the gas and the humidity. 


3. To calculate the integral of the change in enthalpy 
and entropy describing the change in state of air or com- 


bustion gases 


If the state of the gas is known in at least two 
places, the quantities can be calculated which indi- 
cate how the state changes. The most suitable are 


2 To calculate T*—T at given temperature 7 and given 
enthalpy difference Ah, the sub-routine 3.3 described later may, 
for instance, be used. The pressure ratio p*/p can be calculated 
from the entropy difference at constant pressure corresponding 
to these temperature limits, using the sub-routine 3.2, also 
described later. ; 
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eT. <6 (7-7) 
fis = fm - S(p=p) 


i 2 
Ahi+i=5 em (1+ cota) 


= 


|Ahi+ i-Ah;| - € 
BROWN BOVERI! 
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Fig. 3. — Flow diagram for computation of the gas state, 
sub-routine 2 


For explanation, see text. 


the enthalpy difference Ah and the entropy difference 
at constant pressure As,. If the gas may be regarded 
as an ideal gas, these differences can be expressed 
by the following integrals 


et 
As, Lal ie dT 
tT i) D9 SRG & 


Therein c, is the specific heat at constant pressure. 
The enthalpy difference indicates directly by what 
amount the energy of the gas was increased by work 
or the addition of heat, or by a change in the kinetic 
energy. With the aid of the entropy difference the 
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Fig. 4. — Differences in enthalpy and entropy in the enthalpy- 
entropy diagram 
a: For compressors b: For turbines 


For notation, see text. 


polytropic efficiency 


Asr lor, 2” 
R °t he 
1] a ‘. =. 
As, tp aT 
kt TRaear 
Ta 


can, for instance, be calculated (Fig. 4) [2]. The 
subscripts « and w denote that the quantities refer 
to the states « and w. The above definition of the 
efficiency applies to compressors; for turbines the 
reciprocal must be inserted. From this definition 
it will be seen that it is convenient to calculate 
the quotient of entropy difference divided by the 
gas constant, instead of the entropy difference itself. 


326 


THE BROWN BoveERI REVIEW 


VOL. 48, No. 5/6 


BROWN BOVERI 


117217-1 


Fig. 5. — Flow diagram for computation of the specific heat by 
an approximation polynomial function, sub-routine 3.1 


For explanation, see text. 


In addition to these differences, however, it must 
be possible to work out one temperature limit when 
the value of the difference and the other temperature 
limit are known. In many cases an enthalpy differ- 
ence can be determined more accurately from a 
measurement of mechanical power and the mass flow, 
for example, than from a direct measurement of 
the temperature difference. 

The two integrals for the differences in enthalpy 
and entropy depend on the composition of the gas 
and the temperature since the specific heat for a 
particular gas is only a function of the temperature, 
at least within normal pressure limits. When calcu- 
lating by hand, the values of these integrals can be 
taken from the curves applying to a certain gas 
composition or from tables. For the calculation by 
the computer the specific heat at constant pressure 
for a particular gas will be represented by an approxi- 
mation polynomial having the form 


n 


Ra~Hty T+ aT? + ...4,T*= aT 


Closed integration of such a polynomial is possible, 
the values of the two integrals being calculated with 
appropriate constants, as required. 


If the calculation is performed for air on the one 
hand and for the products of a stoichiometric com- 
bustion of definite hydro-carbons on the other, the 
integral values for any mixtures of these gases can 
be linearly interpolated, provided the mixing ratio 
is suitably defined. The calculated values then have 


to be corrected for the moisture contained in the air. 


In the remarks that follow, only one gas will be 
considered, in order to simplify matters. The com- 
puter program must now be able to supply the value 
of the integrals with definite temperature limits, as 
well as one temperature limit for a definite value 
of the integral with a known value for the other 
temperature limit. To cope with these two tasks the 
program is best divided into the following three 


sub-routines. 


3.1 For a given value of the temperature the 
specific heat divided by the gas constant is calculated 
for the dry gas at constant pressure (c,/R)o. This is 
best done by a cyclic process, in which for each run, 
an additional term is included in the polynomial 
having the following form (Fig. 5): 


(4) age RCA n-1) T+ An:2++ «| ine % 
0 


This value finally has to be corrected for the hu- 
midity X. 


3.2 Using a very similar method the values of the 
integrals are calculated for two given temperature 
limits. The choice of the computation procedure for 
the desired integral can, for instance, be governed 
by a quantity z, to which different values must be 
assigned in the two cases. For the calculation of the 
enthalpy integral it must be 0, while for the entropy 
integral it must be | (see Fig. 6). Here too the integ- 
rals have to be corrected for the humidity X. 


3.3 If, for a given enthalpy or entropy difference 
and one given temperature limit, say T,, the other 
limit T,, has to be calculated, a decision must first 
be made (z = 0, z= 1) as to which of the two 
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Fig. 6. — Flow diagram for computation of the integrals of the 
changes in enthalpy and entropy, sub-routine 3.2 


For explanation, see text. 


differences is actually given (Fig. 7). The temperature 
limit 7, is found when the given difference J equals 
the integral F between the given limit 7,, and the 
desired limit 7;,. To find the zero point of the func- 
tion (J — F) dependent upon the desired temper- 
ature it is best to use an iteration. A better approxi- 
mation Tj,; is obtained from the preceding 7; if 
the functional value (J — F) for this temperature is 
divided by the derivative of this function for the 
temperature T;: 


Tigi = T; + 


This is commonly referred to as Newton’s method 
and is ideal in this case because the derivative of 
the integral, i.e. the integrand f, can be calculated 
with sub-routine 3.1, and the integral F with 3.2. 
With appropriate specification of the initial value 7, 
the method converges in the whole temperature 
range. If the improvement in the temperature from 
one step of the iteration to the next is smaller than 
a definite limit ¢, the iteration is discontinued and 
the last calculated approximation for the temper- 
ature limit is the result. 
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Fig. 7. — Computation of one temperature limit when the other 
limit is given, for given enthalpy or entropy difference, sub- 
routine 3.3 


For explanation, see text. 


Storage of Experimentally Determined 
Functions in the Sub-Routines 


Within the calculations described above it is often 
necessary to determine quantities bearing an experi- 
mentally determined relationship to the calculated 
quantities, for instance in sub-routine | the co- 
efficient of discharge «[(d/D)*, D], the expansion 
factor e[ (p,./p,)"/*, (d/D)2] and the density of the 
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water @y (t,); and in sub-routine 3.1 the function 
c,/R (T). It must therefore be possible to determine 
the value of these functions within a certain range 
for any value of the argument. This task frequently 
occurs in other calculations too; examples will there- 
fore be given to demonstrate how it was solved in 
connection with the present sub-routines. 

For non-mechanized calculation such functions are 
usually represented by tables or curves. On changing 
over to machine computation there are two funda- 
mental possibilities for storing such functions within 


the computer program. 


1. Associated pairs of values for the argument and 
the function are stored as a table. In most cases, 
though, the value of the argument for which the 
functional value is desired is somewhere between 
the stored values. Consequently a program must 
also be provided for interpolation of the func- 


tional values. 


2. For each value of the argument the value of the 
function is calculated according to a program, 
provided the relationship can be expressed at least 


approximately as a formula. 


Since the sub-routines are used in very different 
situations, the range in which the functional values 
have to be determined is very large. In order to 
fulfil the requirement for high accuracy of calcula- 
tion and minimum storage capacity under these 
circumstances, the function was stored in every case 
in the form of a program operating on the principle 
of an approximation formula. 

But in every case this approximation formula has 
to be found first. Often it is possible to abstract the 
physical process to be described in such a manner 
that only the significant factors are taken into con- 
sideration, and all insignificant factors disregarded. 
In many cases this physical model can be described 
mathematically. Such formulae expressing the physi- 
cal process direct can then be utilized as the basis 
of the program. 

An example is the calculation of the expansion 
factor for sub-routine 1. According to the VDI 


definition, ¢ is given by the formula 


Panga: 1 — pin? P? x 
do 1 — wim? P? Le P* ied 


(1 — P*4) 


in which 

/4, = contraction factor for gases 

%, = discharge coefficient for compressible gases 

«= discharge coefficient for incompressible liquids 


*=ratio of the specific heat at constant pressure 
to that at constant volume 
PEN. 
m= (5) = ratio of the cross-sections of nozzle 


\ 


and pipe 
I 
P= abbreviation for ({:) % 


1 

A closer examination shows that the value of the 
expansion factor is mainly dependent on the quan- 
tity P and on the cross-sectional ratio m, but only 
to a very small extent on the ratio of the specific 
heats which otherwise appear at various points in 
the formula. Thus, for example, with x between 
1-0 and 1.5, the differences in ¢ remain below 0-13 %,. 
If we were now to expand the third factor in powers 
of (1 — P), the ratio of the specific heats can be 
arbitrarily made = | in the various terms of the 
expansion. The series may be discontinued after the 
third term without the accuracy suffering. Further- 
more, for the standard nozzle, ignoring radial 
expansion, the contraction factor 4, and the ratio 
of the discharge coefficients «,/x% may both be 
equated to 1. Thus the simplified formula 


te 1 — m? 
= V1 — m2 P? 


3 l 2). A eee 
|= Fa-P) +5 0-P P+ ay =P) 


is obtained, which is ideal for cyclic programming. 
The deviation of the values calculated with this for- 
mula from those given by curves plotted from 
measurements and reproduced in the VDI regula- 
tions is less than 0-14°%, provided the cross-sectional 
ratio is restricted to m < 0-6. 

There are, however, other cases in which it is 
either impossible to find a model for the physical 
process, or the formulae derived from a model are 
too complex and cumbersome to permit their use 
in acomputer program. In such cases an attempt must 
be made to find a suitable approximation formula 
which approaches sufficiently closely to the given 
points but bears no direct relationship to the physical 
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process. But there is no point in forcing the function 
to pass through the given points, which are often 
affected by errors of measurement in any case. The 
approximation function would only be rendered 
more complicated by doing so. Using the method 
of the smallest squares it is usually relatively easy 
to find a simple approximation formula. 

This method was used, for instance, to calculate 
the values of the specific heat at constant pressure 
divided by the gas constant ¢c,/R as a function of 
the temperature for sub-routines 3.1. The tabulated 
values (see [3, 4, 5] for example) were mostly 
obtained using statistical formulae of quantum 
mechanics from analysis of the spectra of the various 
gases (e.g. [6 and 7] ). However, these formulae 
are far too complex to be used as a basis for a sub- 
routine in the computation. In contrast, it was 
possible to find a second-order polynomial as an 
approximation to this function for dry air in the 
temperature range 240°K < T < 580 °K, and a 
fourth-order polynomial for dry air and the products 
of a stoichiometric combustion in the range 580 °K < 
T < 1280 °K, the deviations from the tabulated 
values being less than 0-03°% over the entire 
temperature range. Representing the c,/R function 
by approximation formulae has, as already men- 
tioned, the additional advantage of permitting 
closed integration of the polynomials, so that 
the integrals can also be calculated by approxi- 
mation polynomials with appropriate coefficients. 


Conclusions 


The present article has shown how great and 
varied are the advantages gained by the intro- 
duction of digital computation for the evaluation 
of measurements. But it is important not to program 
merely the procedure normally used when calcu- 
lating by hand; rather, the problem must be tackled 
in a way which suits the abilities of the computer. 
This applies particularly to the calculation pro- 
cedure for which, corresponding to the very high 
speed of calculation, cyclic and iterative methods 
are often employed, whereas explicit methods would 
be used for calculation by hand. Likewise for storing 
experimentally determined functions and for the 


build-up of the whole program it is essential to find 
a suitable procedure. To explain this point a des- 
cription of existing sub-routines was given, which 
allow programs to be compiled in a very short time. 
It thus becomes possible to employ a computer 
rationally even for the evaluation of quite small 
series of measurements on different kinds of test 


objects. 
G. D1BEvius 


(KME) E. TARNOCzY 


Summary of the Notation Used in the Text 


Symbol Meaning Units 
A = Cross-sectional area m? 
B= Factor in the formula for mass 

flow kes! m? 
D = Diameter of pipe for measure- 
ment of mass flow m 
F = Integral value [yar eas OF 
Ta dimensionless 
K = Factor in the formula for mass 
flow m2 
M = Mass flow kg/s 
P = Abbreviation for (p,/p,)'/* -— 
R= Gas constant J/kg °K 
T = Absolute temperature ois 
X= Molecular ratio of the water 
vapour in the air to dry air —- 
a= Polynomial coefficient — 
¢, = Specific heat at constant pres- 
sure Jikg “K 
d= Diameter of the standard nozzle m 
f= Function ¢,/R or 3 ¢,|R dimensionless 
or 1/°K 
h= Specific enthalpy J/kg 
k—=Temperature correction for 
nozzle area — 
m= Ratio of areas of nozzle and 
pipe = 
p = Absolute pressure N/m? 
s = Specific entropy J/kg °K 
t= Temperature above the freez- 
ing point of water ra 
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Symbol Meaning Units in = incompressible 
«= Discharge coefficient for stan- ¢ = compressible 
dard nozzle ina pipe roughened p =at constant pressure 
by normal service — 7 =at constant temperature 
é = Expansion factor — W = of water 
: : A = difference 
7 = Polytropic efficiency — 
3 = Pressure detecting element dy- 
namic correction factor — Bibliography 
Se io of ific heats — 
ae ae ESP eae [1] DIN’ 1952, VDI-Durchflussregein, 6th edition, Publ By 
A= Coefficient of superficial ex- VDI-Verlag GmbH., Dusseldorf, Nov. 1948. 
pansion 1/ °C [2] O. ZweireL: The determination of the variation of state in 
(io Contraction factor — turbo-machinery by means of the increase in entropy. 
: Brown Boveri Rey. 1941, Vol. 28, No. 8/9, p. 232-6. 
& = Temperature sensing element 
dynamic correction factor PL [3] Tables of thermal properties of gases. U.S. Department of 
} , Commerce, National Bureau of Standards Circular No. 564, 
o = Density (specific weight) kg/m? 1955. 


Suffixes and symbols 


0 =at ambient temperature 
0 =for dry gas 
a, w=denoting different states 
D = of the nozzle 
1 = before the nozzle 
2 =after the nozzle 


7 = current suffix 


[4] 


[5] 


[6] 


[7] 


C. H. Wu: Tables of thermodynamic properties of combus- 
tion gases. Science Press, Peking 1958. 


J. H. Keenan, J. Kaye: Gas tables. Publ. by John Wiley & 
Sons, New York 1948. 


J. E. Mayer, M. G. Mayer: Statistical mechanics. Publ. 
by John Wiley & Sons, New York 1940. 


H. W. Wot ey et al.: Compilation of thermal properties of 
hydrogen in its various isotopic and ortho-para modifica- 
tions. Journal of Research of the National Bureau of Stan- 
dards, 1948, Vol. 41, RP 1932, p. 379. 
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THE EMPLOYMENT OF DIGITAL COMPUTERS IN THE 
DESIGN OF ELECTRICAL MACHINES 


A useful aid for the calculations involved in the design of 
electrical machines is the digital computer, which can be 
entrusted with frequently repeated routines, when fed with 
appropriate programs. Brown Boveri have been preparing 
such programs since 1957, many of them being repeatedly 


utilized in practice. 


ROGRAMMED ealculation of electrical ma- 
chines is one of the special applications of the 
digital computer. The calculations involved are 
required when designing electrical machines, as well 
as for checking constructional details. Drafts and the 
calculations for such units as transformers, motors or 
generators are continuously recurring in practice 
and take up a great deal of time. A design calcula- 
tion necessitates the employment of a large number 
of formulae and laws, possessing general validity in 
the range of a certain type of design, but which to be 
numerically evaluated individually every time. Owing 
to its ability to perform extensive numerical calcu- 
lations, following a preset program, at high speed, 
the digital computer is an ideal aid for such design 


calculations. 


Design Calculation and its Automation 


As a general rule the design calculation of an 
electrical machine may be divided into three phases: 
dimensioning, determining the properties, and assess- 
ment. When the calculation is carried out by hand, 
these three activities are interlaced to a large extent. 
The skilled design engineer is continually assessing 
the intermediate values and results, with the object 
of recognizing the effect of the dimensions adopted 


681.14:621,.313.001.24 


as quickly as possible. Also the procedure for a 
manual calculation will differ from one case to ano- 
ther. It is Obviously a process which does not follow 
any hard and fast plan. The main object of auto- 
matizing such processes is therefore planning an 
organized general schedule for the calculation. The 
procedure must be moulded into a definite succession 
of rational steps. Furthermore, the task in question, 
if it is to be handed over to a program-controlled 
computer for solution, must be suitable for formu- 
lation. Therefore it is necessary to be able to analyti- 
cally express graphical processes and inevitable esti- 
mates based on experience. 

All these requirements are essential stipulations 
which have to be fulfilled when a design calculation 
by hand is converted into an electronic computation. 
In many cases the task is far from easy, but never- 
theless soluble. Having planned the structure of the 
computation process and formulated the quantities 
to be calculated, the use of the electronic computer 
opens up possibilities offering numerous advantages 


over calculation by hand. 


The Calculation Program 


The program for a design calculation comprises all 
necessary constants and data, as well as the relevant 
instructions, the execution of which in the specified 
order represents the computation procedure for a 
particular type of machine. The mathematical ope- 
rations for which instructions are given process the 
constants of the problem, as contained in the pro- 
gram, and the variable data which are features of the 
design in question. The instructions are so arranged 
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that they can be carried outin the manner appropriate 
to every possible design. The instructions and all data 
not applying exclusively to the particular application 
consequently form a generally valid whole, known as 
the program. Thus each program, since it possesses 
general validity in its field of application, only needs 
to be set up once. The instructions and general 
numerical values are recorded in code form on pun- 
ched tape. The program can be retained in this 
permanent form and is always available for com- 
munication to the stores of the computer. Before 
such a program can be put into action, the variables 
have only to be compiled. 

Depending on the content and the extent, com- 
puter programs may be divided into three categories, 


as follows: 


1. Those which determine the properties of a design 


whose geometric proportions are completely given. 


2. Partly automatic programs which supplement 
designs existing in an incomplete state and deter- 
mine their properties. Partial dimensioning is 
achieved by observing conditions relating to the 
utilization and the behaviour of the finished 
design, and in some cases to prevailing standards, 
and by fulfilling these conditions in the best pos- 


sible manner. 


3. Completely automatic programs, which carry out 
the whole draft in such a manner that the resul- 
tant solution fulfils all stipulations made regarding 


the properties. 


A feature common to all these programs is the 
determination of the properties of the design. The 
difference between the three categories is in the exis- 
tence and the extent of automatic dimensioning. The 
actual process of drawing up the design makes allow- 
ance for, and fulfils as well as it can, the conditions 
which can be stipulated for the design to be supple- 
mented or drawn up from scratch. This process of 
aimed design is also known as optimizing, since it 
aims at providing the best solution, aided by pro- 
grammed assessment of certain results in conjunction 


with appropriate modification of the parameters. 


In accordance with the graduated capabilities, the 
types of programs listed naturally also differ in the 
time taken to set them up, in their structure and 
extent, and finally in their behaviour. With the 
degree of automation, the elegance with which each 
program can be executed increases. Whereas a com- 
pletely automatic program determines the best solu- 
tion in a single sequence, the simple ‘check’ program 
and partly automatic programs have to work out 
several variants before a good solution is found. 
However, the higher the degree of automation, the 
more work setting up involves, the greater the need 
for storage capacity and the greater the volume of 
actual calculation. The flexibility of a computer 
program is influenced in the opposite sense. With 
programs having a lower degree of automation, the 
flexibility is greater, because with them the suitabi- 
lity of the solution is left to human judgement. With- 
out changing the program it is possible for the manner 
in which the assessment is made to change. The 
chance of changing a program also diminishes with 
increasing automation. Completely automatic pro- 
grams are consequently only employed for revo- 


lutionary types of machines. 


Preparation and Employment of a : 
Design Program 


In Fig. 1 the steps leading up the completion of a 
program, and during its operational employment are 


depicted schematically. The upper leg in this diagram 


LE aE 


refers to the preparation of the program. In a single 

operation the calculation formulae and all necessary 
mathematical and technical constants are translated ~ 
into the language of the computer. The now encoded | 


program is next punched on the tape with the aid 


of a sheet printer which also acts as a punch. Thus, at 
the same time as the punched tape is being produced, 
a record in clear is typed out showing the contents of 
the tape. Only when the tape is finished is it possible 
to impart the instructions and general data to the 
stores of the computer. Since the program possesses 


general validity within its field of application, it is 
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Fig. 1. — Diagram showing the preparation and employment of a design program for the digital computer 


The upper leg shows the preparation of the general computer program, which is stored in the form of punched tapes in the library 


of the Computer Centre. The lower part of the diagram shows how the program is fed into the computer and the problem calculated. 


The design engineer collects the various performance data of the machine to be calculated, which are read into the computer to- 


gether with the existing program, all in the form of punched tape. On completion of the computation, the results are passed back 
to the design department. 


kept in the library of the Computer Centre in the 
form of a punched tape and is thus available for use 
again at any time. 

The lower leg of the diagram in Fig. | applies to 
the preparation and execution of the calculation. 
The responsible engineer from the design department 
compiles the variable data for one or more design 
variants. These data contain the ratings of the electri- 
cal machine being calculated and, depending on the 
nature of the program, further requirements relating 
to the form of the design or the attainment of certain 
properties. The variable data are transferred to a 
punched tape by a sheet printer, and that completes 
the preparations for the calculation. 

Only the succeeding operations are performed in 
the Computer Centre. They can be performed by 
the operator without any assistance by the responsible 
engineer. The program kept in the library and the 
variable data are fed into the computer store 
through input elements. Then the computer performs 
the instructed operations with the given numerical 
values, the result appearing as a punched tape. 
Independent of the computer the results can be 
translated into clear by the sheet printer. In this 
form they can be handed over to the design de- 
partment. 


The Uses and 


Limitations of Electronic Computation 


Program-controlled design calculation offers a 
variety of advantages over the classical method of 
calculating by hand. When setting up a computer 
program the methods of manual calculation are not 
exclusively translated into the computer language, as 
a general rule. In view of the high speed and the large 
storage capacity, the amount of actual calculation 
work can be greatly increased. The programmed task 
may therefore be based on more extensive and refined 
methods of solution, which benefits the calculation, 
in that it is thereby made more comprehensive and 
more accurate. In many cases the use of the digital 
computer has made possible the employment of 
methods involving far too much work to allow them 
to be used in manual calculations. 

By means of the program a great deal of mathema- 
tical work is transferred to the computer, which 
previously was the entire responsibility of the design 
engineer. This work is performed rapidly and free 
from errors. By way of example, the normal calcula- 
tion of a turbo-alternator used to occupy two full 
working days, when done by hand. With one of our 


semi-automatic programs this work can be completed 
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within five minutes, disregarding the time taken to 
feed in the program. The engineer, now largely 
relieved of the burden of exhaustive calculations, is 
thus able to devote more time to a thorough criticism 
of the results. To judge by experience gained so far, 
he utilizes the greater part of time saved in thinking 
out refinements to the design, so that changing over 
to electronic calculation does not represent any great 
reduction in the number of design engineers, but an 
improvement in the quality of the calculated product. 
In fact, the high speed of the digital computer allows 
several variants to be calculated in a very short time. 
It provides a means of selecting the best design from 
several in every case; this was not always possible 
with manual calculation owing to the time and the 
amount of work involved. With transformers and 
rotating electrical machines it has been established 
that electronic computation has resulted in a appre- 
ciable improvement in the quality. This may be 
expressed either by an improvement in the opera- 
ting properties or a reduction in the weight, depen- 
ding on which is considered more valuable. 

Generally speaking, it may be claimed that the 
change from manual to electronic calculation has 
opened up numerous new possibilities. Moreover, 
design programs are very versatile in their appli- 
cation, and suitable for work in connection with 
tenders, isolated new designs and the compilation of 
price lists. 

A technical limit to the scope of electronic design 
calculations is the storage capacity of the particular 
computer. Now one of the features of design pro- 
grams is the large number of instructions and logical 
operations, and the variety of the data required and 
individual results. The storage capacity needed for 
our programs, as set up for the Siemens 2002 com- 
puter for instance, is between 3000 and 11 000 cells, 
the largest capacity being required for the comple- 
tely automatic programs. Thus, to undertake design 
calculations without any limitations on the scope of 
the task it is essential to employ digital computers 


having a large storage capacity. 


Economics of Electronic Computation 


In this connection it is equally important to con- 
sider the economics of electronic design calculation. 


The work involved in setting up a program-controlled 


calculation must be compared with the attainable 
usefulness. It must be pointed out that electronic 
computation, with its associated advantages, is only 
made possible by the fairly expensive business of 
preparing the program. Depending on the degree of 
automation, it may take from 500 to 3000 man- 
hours — a not inconsiderable figure — to prepare a 
single program, quite a high proportion of which can 
be left to suitably trained auxiliary personnel though. 
The fixed costs must be offset against the attainable 
improvement in the quality of the product and 
reduction of its price, in addition to the saving 
effected in the engineer’s time, if electronic compu- 
tation is to justify itself as an economical process. 
But the more often a program is utilized, the smaller 
each variant’s share of the initial costs becomes. An 
economic gain is therefore assured when the extent 
of a program is adapted to the frequency with which 
it will be used and the possible improvement in the 
quality of the calculated product. In contrast, the 
use of an electronic computer is uneconomical when 
the cost of preparing a program exceeds the cost of 


calculation with simpler means. 


Development and Present State 
of Electronic Design Computation at 
Brown Boveri 


Commencing in about 1957, exhaustive studies 
were carried out to investigate the possibilities of 
using program-controlled computers for calculations 
needed for the design of electrical machines, a num- 
ber of different computers being tried out to establish 
their suitability. In the course of these investigations 
and subsequently our first design programs were 
prepared for the different types of computer which 
were available. Even at that time some of these 
programs were usable for practical calculations. 
The experience gained then proved most useful 
when it came to the preparation of programs for our 
own computer. An early start was made with the 
conversion of existing programs to suit the computer 
decided upon, and to prepare new programs for it. 
Thus by the time the digital computer was ready to 
commence service several design programs were 
already available for practical utilization. Since then, 
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of course, others have been added. Now the library 
contains tried programs in different degrees of auto- 
mation for numerous electrical machine designs. 
They are constantly being used for practical pur- 
poses and have rendered valuable service in the 
preparation of tenders, isolated new designs and the 
compilation of price lists. 


Induction machines 


In the field of induction machines partly automatic 
and check programs are being employed with great 
success. Programs for slip-ring machines, squirrel- 
cage motors with double-cage rotor, or with deep- 
slot rotors, have existed for some time. In these 
programs, apart from determining the electrical, 
magnetic and thermal conditions, stress is laid on the 
determination of the operational performance. A 
program serving this purpose will be described 
briefly in the next chapter. 


Synchronous machines 


For the field of synchronous machines there are 
efficient, partly automatic programs for salient-pole 
generators of low output and for large hydrogen- 
cooled turbo-alternators. Apart from their use for 
the automatic dimensioning, these programs are 
employed to calculate all the major properties. For 
the high-power machines in particular, a detailed 
price and weight calculation is of great economic 
value. It has therefore been included in the pro- 
gram for turbo-alternators and permits a comprehen- 


sive assessment to be made of the solution. 


Transformers 


For small transformers with ratings of up to 
650 kVA a program has been drawn up with the 
main object of gaining experience. It can be used for 
price list calculations and for variants on a basic 
design. It was used in the development of a new 
range of units. Afterwards a new program was 
devised to cover the whole range from 20 to6000kVA, 
provided no regulating winding is involved. It is at 
present being introduced, both for the development 
of new series and for dealing with quotations and 
orders. 


q 


Apart from these typical design programs, there 
are also several programs for special purposes, such 
as the calculation of short-circuit forces,! magnetic 
fields, the current distribution in parallel groups of 
coils, and the distribution of impulse voltages. 

At the present moment the program-controlled 
calculation of electrical machines is being rapidly 
expanded at Brown Boveri. Numerous other pro- 
grams with different degrees of automation are being 
prepared. From the success gained so far it is evident 
that electronic computation will in due course find 
its way into all branches of engineering design within 


the Company. 


Description of a Partly Automatic Program 
for Calculating Squirrel-Cage Motors 


Task of the Program 


With the program it should be possible to deter- 
mine the electrical, magnetic, thermal and starting 
properties of any squirrel-cage motor of given 
constructional design. In the event of the design 
department not prescribing a definite stator wind- 
ing, the program should also be able to find the best 


stator winding and carry out the calculation for it. 


Principle of the Program 


When setting up a program, a block diagram is 
first prepared, showing the various stages of the 
calculation in very general terms. Partial problems 
which are complete in themselves are denoted by a 
box containing a keyword (see Fig. 2). The next 
step consists of working out these partial problems 
in detail, thereby producing a flow diagram. This 
is then a complete set of instructions which finally has 
only to be translated into the language of the elec- 
tronic computer. Let us now examine the block 
diagram, with reference to Fig. 2: On starting the 
computation the motor data written on punched 
tape are first transferred to the store of the computer. 
They contain the stipulated operating conditions 
and the constructional data of the motor. The first 


1M. CurisrorreL: Calculation of short-circuit stresses in 
transformer windings with the aid of digital computers. Brown 
Boveri Rev. 1960, Vol. 47, No. 5/6, p. 321-8. 
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Read in motor data 


Yes 


bas ea; Is the stator winding prescribed ? rigs cage Design a stator winding 


Calculate the leakage due to harmonics 


for stator and rotor 


Estimate the leakage factor 7 


Estimate the factor 7 
for deformation of the magnetic field 


Calculate the magnetic field 
in the iron 
Calculate 7 
Was 7 estimated 
of motor accurately enough? 


Was T estimated : : 
B Estimate the efficiency 7 
accurately enough? 


Calculate the losses and Heyland 


diagram 


Calculate temp. Yes Was 7 estimated A 
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Galculate price Design a new stator 

: : winding 
Calculate starting 

characteristic Print out results 

y * = A 3 Yes 
Was the stator winding specified ? Recall stored results 


Does this winding produce the sii atiafor ill wore drone 


best possible motor? 
Fig. 2. — Block diagram for the computation of squirrel-cage motors 


Calculate leakage 


Have all possible stator 
windings been calculated } No 


117276/le 


BROWN BOVERI 


. 


Rectangular boxes denote closed sub-routines. Those with rounded ends denote decisions to be made by the computer. Note the 
links at three points. These indicate that parts of the calculation have to be repeated several times because a certain factor is unknown 
to begin with and must consequently be estimated. Finally, when the factor can be calculated, the accuracy of the estimate can be 
judged. If the accuracy is inadequate, an improved value is estimated with the aid of the calculated factor and the calculation 


repeated. 
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step is to decide whether these data determine the 
stator winding. If they do, the design of the motor is 
already fixed and the task of the program is solely to 
determine its properties. If not, the computer has to 
decide on the data for the missing winding, according 
to a separate instruction, and only then can it carry 
on with the determination of the properties. First a 
Fourier series has to be evaluated to obtain the 
leakage due to harmonics of the stator, followed by 
another infinite series to determine the same leakage 


for the rotor. 


To enable the magnetic field in the motor to be 
carefully investigated, and from it the no-load and 
short-circuit currents, two coefficients y and Tt are 
required, whose magnitude depends on the satura- 
tion of the iron and the constructional design. The 
stray factor tT summarily allows for the stray field 
in the motor; the factor y describes the deviation of 
the field from the sinusoidal form when the iron is 
heavily saturated. Neither of them is known to begin 


with and they are both obtained in the course of the . 


calculation. They must therefore be estimated to 
start with. Having calculated the fild in the iron, y 
can be worked out and compared with the estimated 
value. If the deviation is very great, a better value is 
assumed for y and the calculation repeated. This 
procedure is repeated until the estimate and the cal- 
culated value agree sufficiently closely. Not until 
then can the stray field, and from it the factor tT, be 
calculated. The same procedure is adopted as for y, 
and for the efficiency 7. Finally the starting charac- 
teristics, i.e. the curves of stator current and torque 
against the speed, are computed. If the stator winding 
was given, the calculation is now ended and the 
results are simply printed out in a list. 

However, if the designer is free to use a stator 
winding of his own preference, he can leave the choice 
of this to the computer by not prescribing any rele- 
vant data in the program. The program then ensures 
that all possible forms of winding are considered 
successively, and the properties of each motor cal- 
culated. From the finite number of such solutions, the 
computer then has to select the best by comparing 
the price and the performance figures, such as tempe- 
rature rise, efficiency, power factor, starting and 
stalling torque. Every time a better solution than all 
preceding ones has been found, the results are trans- 


ferred to a place specially provided in the store. Here 
is always kept the best solution determined and, after 
calculating all possibilities, the data for the best 
motor. 

If guarantee values are prescribed, these and the 
motor data are fed into the computer, which first 
checks whether they can be fulfilled. If any of the 
guarantee stipulations is not met, the results of the 
calculation are not stored and are erased on commen- 
cement of the next calculation. It is, of course, 
possible for the conditions stipulated for guarantees 
to be so severe that there is no winding which can 
possibly fulfil them. In this case the computer prints 
out a short notice instead of results, indicating to the 
engineer all the information he needs to pass jud- 
gement on his design. 

When calculating motors by hand, graphical 
methods (such as the Heyland diagram) are used 
which, with the aid of analytical geometry, can be 
made readily accessible for the electronic computer. 
The design engineer also employs numerous curves 
and tables which can also be provided in an analytical 


form by various means. 


Practical Work with the Program 


Having set up the program and checked it, the 
engineer deposits it in the Computer Centre with 
his instructions for its use. Each design engineer, in 
addition to receiving the set of instructions, also has 
a flow diagram showing the method adopted for 
computation, and a set of printed lists for input data, 
which he has to fill out each time before utilizing the 
program and hand over to the Computer Centre. 

On receipt of an order or enquiry the design engi- 
neer first prepares a draft that can be easily calculated 
from existing types, from motors which have already 
been built, and utilizing his own experience. The 
data of his draft are entered in on a printed form, 
handed over to the Computer Centre, where they 
are transferred to punched tape and subsequently 
computed. All the work in the Centre is carried out 
by auxiliary personnel. It is of course possible for 
the engineer to remain in the Centre and observe the 
results as they come out of the computer, permitting 
him to repeat the calculation with different data. 
Without being present in the Centre he can also 


) 
oO 
© 


leave instructions for several sets of data to be calcula- 
ted and can compare the results afterwards. In this 
manner the design engineer can draft motors quickly 
whose details are matched in every respect, a process 
which would be quite impossible with the slide-rule, 
above all for very small machines. He is thus relieved 
of almost all manual work, and can concentrate his 
efforts on the utilization of his experience. This com- 
bined operation by the design engineer and the com- 
puter has proved very successful. 

The time taken to feed in the design program des- 
cribed is in the region of six minutes. The calculation 


of one kind of winding lasts about three seconds. If 
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the computer has to investigate one hundred different 
windings for a particular motor it requires roughly 
five minutes for this. Printing out the complete list of 
results, involving some 200-300 numbers and con- 
taining the prescribed motor data, lasts barely two 
minutes. A skilled engineer requires approximately 
one day for such a calculation. The program descri- 
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bed is therefore an ideal means of working out — 


several different designs of an induction motor and, — 


by comparing the results obtained from each calcu-_ 


lation, to find the best solution. 


G. NEIDHOFER 
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IMPROVING THE SMOOTH-RUNNING QUALITIES 
OF TURBO-ALTERNATORS 


The author describes the employment of the digital com- 
puter during the design of turbo-alternator windings con- 
sisting of bars made of hollow-drawn conductors. Irregulari- 
ties arising during the manufacture of the conductors are 
balanced by using a computed bar assembly, thus eliminating 
one of the main thermal factors affecting the smooth-running 


qualities of such rotors. 


Temperature Rise and Smooth-Running 
Qualities 


UITE INSIGNIFICANT differences 
temperature of diametrically opposed zones 


in the 


in the iron of turbo-alternator rotors can produce 
distortion seriously affecting the smooth running of 
the machine. Such deviations from the normal tem- 
perature rise in operation are éaused by thermal dis- 
turbances originating in the excitation winding situ- 
ated in the rotor slots. The generation of unequal 
amounts of heat in the slot filling due to the passage 
of current, and the different rates of dissipation of 
the heat from the copper are the main disturbances. 
When electrically loaded, the local temperature rise 
of the rotor iron is always influenced by the adjacent 
slot. If the passage of current through diametrically 
opposed groups of slots produces different temper- 
ature rises in the surrounding iron, this results in dis- 
tortion of the rotor body which, when running at 
speed, is apparent as a deterioration in the smooth 


running of the alternator. 


The Hollow Conductor 
as a Source of Disturbance 


In the Company’s turbo-alternators for medium 
and high outputs most of the active part of the rotor 
winding consists of hollow copper conductors, the 
coolant flowing through the cavity and the metal 


681.14:621.313.322-81 


acting as the conductor of electricity. With this direct 
system of cooling the rotor iron is in thermal contact 
with the insulating lining, and thus indirectly with 
the bars in which the heat is generated when the 
machine is running. Any differences between the 
hollow conductors allocated to the various slots can 
therefore lead to unequal peripheral temperatures 
in the iron when current is flowing through the 
conductor. 

Owing to the irregularities inherent in the manu- 
facturing process, there are indeed differences be- 
tween individual hollow-drawn copper sections, 
which may adversely affect the temperature rise in 
service. On account of the essential manufacturing 
tolerance, the cross-section of bars composed of such 
conductors (hereafter referred to as hollow bars) may 
consequently not be exactly alike, despite theoretical 
equality of their dimensions. Whereas the external 
dimensions of the bars vary in very close limits, 
the ratio of copper cross-section to cavity is not 
always constant. Differences between the actual 
contour of the cavity from one bar to another are 
quite small, and the roughness of the surface of 
the duct is practically the same in all bars. Along the 
length of a single conductor the cross-section varies 
imperceptibly, so that the differences originating 
during manufacture can only be expressed by the 
ratio of the duct area to the copper cross-section. 

An improvement in the dimensional accuracy of 
the drawing process, which is already difficult enough, 
could only be achieved by increasing the cost out of 
all proportion. The measures described later indicate 
that such a step is not even necessary, as it has proved 
possible to counteract the irregularities mentioned 
above. 

According to the foregoing remarks and owing to 
the strict dimensional accuracy of the exterior of the 
hollow conductors, it is sufficient to regard either the 
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Fig. 1. — Temperature rise of a hollow conductor carrying 
current and cooled by a stream of gas, shown in terms of its 
cavity area Ax and the copper cross-section Acu 


Zone I: 
Zone II: 


low copper losses, heat dissipation too small 
moderate to high copper losses, adequate heat dissi- 
pation 


Zone III: very high copper losses, strong but still inadequate 


heat dissipation 


area of the gas duct or the copper cross-section as the 
sole variable property of the bar governed by its 
manufacture. Since the internal cross-sectional area 
of a tube can be determined directly with ease, and 
with adequate accuracy, every hollow conductor is 
hereafter characterized by the cross-sectional area 


of its cavity. 


Cavity Cross-Section as Characteristic of the Bar 


Differences in the size of the copper cross-section 
and the area of the gas duct affect the heat generated 
by the passage of current and the rate of cooling by 
the gas circuit. Two components—the rate at which 
the heat arrives at the surface and the rate of flow 
of the coolant through the duct—are together respon- 
sible for the temperature rise in a hollow conductor, 
but for such studies it may, under certain circum- 
stances, be sufficient to know only one of the two 
cross-sections. There is in fact an interdependent 
relationship between the temperature rise in a hollow 
conductor and the size of its cavity, when the external 


dimensions remain unchanged and all other quan- 


tities in the system are constant. This relationship is 
clearly illustrated in Fig. 1. It is important to estab- 
lish that, depending on the size of the cavity area, 
any change in it can cause either a positive or nega- 
tive change in the temperature rise. Having regard 
to the excitation losses which are involved, the con- 
ductors of Brown Boveri alternators are designed 
to be within the left-hand section of the curve and 
not too close to the temperature minimum, so that 
for all applications the temperature change is in the 
same sense for all changes in the cavity area. On 
this assumption, and assuming there is not an undue 
amount of scatter in the cavity area from one tube 
to another, it is permissible to assume a linear 
relationship between changes in the cavity area and 
changes in the temperature rise, the relevant section 
of the curve being replaced by its tangent (see Fig. 1). 
Hence, under the above assumptions, it is admissible 
to dispense with the calculation of the temperature 


rise of every one of the large number of different 


hollow conductors, and to evaluate the bars alone | 


on the basis of the deviation of their cavity area from 


an existing mean value. Furthermore, it is possible 


and sufficiently accurate to measure the quality of — 


a slot filling according to its mean cavity area devi- 


ation because, according to the design on the draw-_ 


ing, every bar in the slot is supposed to possess almost _ 


the same thermal weighting. 


Preventive Measures 


It is possible to distribute the hollow conductors - 


of different cross-section in such a manner between 


the individual slots that their different thermal effects 


cancel each other out, as nearly as possible. A pre- 


requisite condition for this is the measurement of the — 


characteristic value for every bar to be incorporated. | 


By means of a systematic process of combination, — 


groups of bars are formed which, according to the 
measurements obtained, should exhibit nearly the 
same behaviour when heated. The thermal effect 


of the residual deviations between the groups so 


formed can be almost eliminated by placing groups j 


of bars having almost the same properties in dia- 
metrically opposed slots. 

The method as such represents a precautionary 
measure which almost completely eliminates the 
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effect of disturbances resulting from manufacture, 


thereby rendering a major contribution to the 


improvement of the smooth-running properties of 


rotors equipped with hollow bars. By suitably organiz- 
ing the process, its execution in practice involves very 


little extra time or work. 


Description of the Calculation Process 


The mathematical processing of the measurements 
obtained for the bars and carrying out the process 
of combination referred to above, together form a 
task, for the solution of which the digital computer 
is ideal. At Brown Boveri a routine has been used for 
some time, which was specially devised for the alloca- 
tion of bars to the rotors of turbo-alternators in the 
course of manufacture. The computer is fed with the 
individual measurements and, apart from providing 
the quality factor for each conductor, yields a com- 
plete insertion programme which can be used by 
the workshop. 


Principle of the Computer Program 


Fig. 2 shows a block diagram illustrating the 
overall structure and the principle of the computer 
program. The remarks which follow are intended 
to give an explanation of certain parts of the prob- 
lem—amarked with numbers in boxes in the diagram. 
To demonstrate the process some of the more im- 
portant steps are illustrated in Fig. 3. The term 
“item” used in the block diagram refers to all hollow 
conductors of the same cross-section and the same 
length. The slot filling of the turbo-alternators in fact 


comprises various conductors of different sizes. 


Computation of bars (1)! 


In box | the stored results of the measurements 
of the hollow conductors are evaluated and all neces- 
Sary preparations made for the subsequent process 
of allocation of the bars. 

The calculation of the characteristic value for each 
bar—in this case the cavity area—depends on the 
method used for measurement. Having regard to 


_1 The bracketed numbers refer to the similarly numbered 


= in Fig. 2. 
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later operations, the bars are arranged according to 
the size of their cavity areas, and then stored. If the 
number of bars delivered exceeds the number re- 
quired for assembly, first those bars are sorted out 
whose absolute cavity area deviation from the exist- 
ing mean is extremely large. Next the mean area of 
the bars to be used is calculated. The above calcula- 


tions are performed for each item (see also Fig. 3a). 


First formation of groups of bars (2) 


From the succession of bars, as they are now ar- 
ranged, groups are formed which are later allocated 


to fill the respective slots (see also Fig. 3b). 


In order, to attain a good combination of bars 
without spending too long on the calculation, it is 
essential not to leave the first allocation of the bars to 
chance, but to carry it out according to a definite 
plan. The quality of each group is determined by 
the relative deviation of its cavity area from the mean 


value calculated for all the bars used. 


Improving the combination (3) 


The object of this is to change the momentary 
combinations in their composition, in such a way 
that the deviations in their area are decreasing in 


magnitude. This is done as follows: 


From the group of bars whose deviation has the 
largest absolute value, which may be referred to as 
the exchange group, the bar occupying a particular 
layer is tried in the corresponding layer in other 
groups (Fig. 3b). For each such experiment new 
deviations are obtained for both groups involved 
in the exchange, if the cavity areas of the inter- 
changed bars differ from one another. Only those 
experiments are kept on record in which the absolute 
value of the new deviations are smaller than that of 
the extreme exchange group before the experiment. 
From the remaining opportunities for exchange each 
time, the one giving the best quality will be selected. 
This process is repeated for each layer until no 


further satisfactory exchange can be made. 


Allocation of groups of bars to rotor slots (4) 


All that remains now is to determine where the 
groups of bars visualized as slot fillings shall be 
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Read in measurements 


First bar position 


Calculate cavity areas 


Arrange numbered bars in order of cavity area 
Eliminate superfluous bars 


Actual mean cavity area 


Output of number and cavity area of each bar with its 


deviation from desired and mean values 


All items dealt with 


Yes 


First formation of groups 


Average area deviation of each group 


Extreme group as exchange group 


Sort out best exchange with one of the 
other groups 


Retain result, if exchange good 


No 
Yes 


Arrangement of groups according to average 
deviation of cavity area 


Allocation of groups to rotor slots 


Print out assembly instructions for the groups 


of bars, showing deviation of each slot filling 


BROWN BOVERI @ 


incorporated. The object of guided insertion of the 


bar groups in the slots (Fig. 3c) is to eliminate, as far 
as possible, the effect of the small deviations between 
the groups which still remain in spite of the im- 
provement process. This is achieved by placing 
groups of bars having almost equal quality factors 


in diametrically opposite slots. 


Next item 


116366 “Ne 


Fig. 2. — Block diagram of the com- 
puter program ; 


Issuing the assembly instructions (5) | 


The computer prints out in clear the calculated — 
composition of the slot fillings, their respective devi- — 
ations and the allocated slot number. The resultant 


chart is then sent direct to the workshop as assembly 
instructions, according to which the assembly of the 


rotor winding can proceed. 
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Fig. 3. — Diagram illustrating the principle of aimed bar 
distribution 


a: Individual bars arranged according to item number and 
deviation 
b: Groups of bars and lines of exchange during the improvement 
process in one layer 
c: Allocation of groups of bars to rotor slots 
The number of bars per slot and the number of bar items have 
been reduced for the illustration. 


Effectiveness of the Program 


The efficiency of the combination process de- 
scribed in the foregoing can easily be judged by com- 
paring the stages before the formation of the groups 
and after they have been sorted out. Fig. 4 shows two 
distribution curves from a case which actually oc- 
curred in practice, from which the initial state and 
the effect of planned bar combination can be com- 
pared. Curve a shows the distribution of the indi- 
vidual bars as a function of their deviations, as 
measured. From a total of 576 bars it is required to 


- 
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Fig. 4. — Distribution curves from an example in practice 


h = Relative frequency 

€p = Relative deviation of cavity area 

Curve a: Distribution curve of the individual conductors 
Curve b: Distribution curve of the slot fillings 


form 64 groups. The r.m.s. error of the groups 
formed in the first step (Fig. 2, box 2) was about one- 
fifth of that of the individual bars. During the subse- 
quent process of improvement (Fig. 2, box 3) the 
mean error in the cavity area was again reduced by 
96 °,. The distribution curve of the improved groups 
is shown as curve b in Fig. 4. Its mean error is 
extremely small. The digital computer goes one step 
further by arranging the groups of bars in the rotor 
slots in such a manner that the slot fillings cannot 
exert any further thermal disturbance effect. 

A feature of the computer program described is 
that it contains rather more logical operations that 
arithmetical. The program comprises roughly 2600 
instructions and, when processing the maximum 
permitted amount of information, requires about 
5600 storage cells. Depending on the number of bars 
and the quality of the conductors, the time taken to 
complete the calculation may be between 70 and 
120 minutes, which includes the process of printing 
out the results. 

A suitable organization has been devised for the 
employment of the program and for punching 
the various measurements, with the result that the 
winding shop can be given the mathematically 
computed assembly instructions at the latest one day 
after receipt of the measurements of the hollow 
conductors. In this case the digital computer is 
rendering invaluable service in the course of a 
manufacturing process. 
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THE CALCULATION OF ELECTROMAGNETIC TRANSIENTS 
ON LINES BY MEANS OF A DIGITAL COMPUTER 


The calculation of discontinuous transients on electric 
transmission lines and cables is very laborious in most cases. 
Following some introductory remarks on typical difficulties, 
the authors explain how the Schnyder-Bergeron method, 
originally devised as a graphical method, can be modified 
to suit the requirements of digital computers. Various 
examples which have been worked out with an appropriate 


program confirm the practical applicability of the method. 


1. The High-Frequency Transient 


N ORDER that all machines and apparatus con- 

cerned with the transmission of electricity may 
be safely and rationally dimensioned, it is most 
important to know something about the over- 
voltages likely to be encountered. Herein a special 
role is played by the high-frequency overvoltages 
which occur after lightning strokes and certain 
switching operations. A characteristic feature of such 
phenomena is that the electrodynamic disturbance 
of the equilibrium of the lines appears in a wave form. 

If full account is to be taken of this special charac- 
ter of the phenomena, the lines and cables must be 
regarded as a continuous structure with distributed 
constants. This applies when either analogue or 
digital computers are used. The strict mathematical 
treatment of continuous structures, such as lines, 
is known to be rather difficult, and almost impossible 
on an analogue computer. 

By substituting chains of quadripoles for the lines 
in the calculation it is possible to find a way round 
the difficulties which prevent them from being 
regarded as continuous structures. In particular, the 
calculation on an analogue basis with the aid of the 
network analyser becomes feasible. But this pro- 
cedure suffers from the serious drawback that, when 


travelling waves with a very steep wave-front have 


681.14: 621.315.064 


to be investigated, there is a risk of the calculation — 


or the analogue simulation of the conditions at 
the wave-front being wrongly reproduced. 

The method described below is based on the 
exact equations of the line with no losses and in 
simple cases is ideal for programming for the digital 


computer. 


2. Various Possible Mathematical Methods 
of Solution 


The equation of the line with no losses is the 
simplest partial differential equation of the second 
order of the hyperbolic type, indeed it is the first 
partial differential equation ever used [1]. D’Alem- 
bert (1717-83) and Daniel Bernoulli (1700-82) found 
two quite different ways of solving this equation, 
which is simultaneously the equation of the vibrating 
string. Bernoulli’s solutions are founded on the eigen 
function and are comparable with the Fourier series. 
In contrast D’Alembert’s method utilizes the wave 
character of the solutions of the hyperbolic differ- 
ential equations. 

It is one of the peculiarities of the hyperbolic 
differential equations that they can assume dis- 
continuous solutions corresponding physically to 


travelling waves with an infinitely steep wave-front. — 


The expansion of solutions according to Fourier 
series is not very suitable for representation of dis- 


continuous functions because the series do not con- 


verge easily, and at the points of discontinuity do — 


not properly reproduce the solutions, despite the 
fact that these solutions can formally be regarded 
as exact (Gibb’s phenomenon) [2]. Moreover, a 
similar phenomenon is observed when, for the cal- 
culation, the line is replaced by a chain of quadripoles. 
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An exact reproduction of travelling waves with dis- 
continuous points cannot be attained, even with 
very close sub-division. D’Alembert’s method avoids 
this difficulty and is used as the basis of the method 
adopted herein. 

For hydraulics the 
method was first devised by Schnyder [3] and 


Bergeron [4], working quite independently, and was 


calculations in technical 


visualized as a graphical method for the calculation 
of transients in penstocks. Programs for computation 
by digital computers already exist for such applica- 
tions [5].! 

These methods, devised by engineers, adopt a 
rather odd system of mathematical expression, but 
it is quite clear and will therefore be retained in this 
article. From the mathematical aspect it may be 
regarded as an application of the theory of character- 
istics, and it is possible to derive the basic formulae 
used from Riemann’s method of integration for 
hyperbolic differential equations. 


3. Fundamentals 
of the Schnyder-Bergeron Method 


3.1 The Line Equations and D’ Alembert’s Solution 


The following notation is used: 
t= time 
x= co-ordinate on the line 


7=current, positive in the direction of positive co- 


ordinates 
“= voltage to earth 
L=inductance per unit length 
C= capacitance per unit length 


z=|/7 = impedance of the line 


l 
Be yic = speed of propagation of the pheno- 


mena along the line 


For the line with no losses the following equations 
may be written 


1 The method given in the quoted paper is, however, based 
on the theory put forward by Allievi, which also relies on pro- 
gressing waves, but which is purely analytical. 


_ Ou ons 00 
Ox ol (1) 
Sf C Ou 
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By eliminating one of the two dependent variables, 


we obtain: 
OF2 1S O22 
Ox? c* Ot? & (2) 
Oru Al 8G ee 0 
Ox? an Le 


F (») and f(y) are any two functions of the variable y. 


‘The two functions 


u(x,t) = F(w—ct)+ f(x + ct) 
u (x,t) = ZF (x — ct) — Zf (x +4 ct) 


(3a) 
(3b) 


separately satisfy the appropriate equation from (2) 
and both satisfy the system (1); they yield the general 
solution to the problem. F(x — ct) describes a wave 
advancing in the positive x direction while f(« +- ct) 
describes one travelling in the negative x direction? 


(D’Alembert’s solution). 


3.2 The Fundamentals of the Schnyder-Bergeron Theory 


Equation (3a) is multiplied throughout by Z, and 
then successively added to and subtracted from (3b), 
giving: 

u(x,t) + Zi (x, t) = 
u(x,t) — Zi (x, t) 


2ZF (x — ct) 
—2Zf (x + ct) 


(4a) 
(4b) 


| 


It may be observed that, in equation (4a), when 
the left-hand side 
constant; the same applies to the left-hand side of 
(4b) when x + ct is constant.3 


For an observer moving in the positive direction 


x — ct is constant, becomes 


with the velocity c, the expression u(x, /) + Zzi(x, t) 
accordingly appears constant. The same applies to 
the expression u(x,t) — Zi(x, t) for an observer 
in the opposite direction. 

In an i, u plane (Fig. 1) wu + Zi = constant corre- 


sponds to the straight line g, and wu — Zi = con- 


2 The theory of hyperbolic differential equations shows that 
F and f may also be discontinuous [6, 7]. 
3x + ct =const. and x — ct = const. are the character- 


istics of the differential equation [6, 7]. 
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stant to the line g_ with the same slope in the opposite 
sense, the sign referring to the direction in which 


the observer is travelling. 


On the line n+ 1 points are given, namely 
Vas esi APE 


(Fig. 2). x) = 0 and x, = / (the length of the line) 


coincide with the two end points Py and P,,. It is also 


..P,, having the co-ordinates X%9...X;.. .X_ 


assumed that 


My << .-- ne 


t; denotes the transit time for the wave over the 
distance P;_, P;. 


(5) 


— x4) /¢ 


Let us now consider three successive points P;_,, 
P; and P;,, on the line, assuming provisionally that 
P; is an internal point). 
In this case it is presumed that all t; are whole- 


j>2i1 andj <n—-1 (ie. 


number multiples of the basic time Tp): 


Tj — m; To 


In the remarks which follow, only whole-number 


multiples of tT) will be considered 
t= MT, 
using m as a new time-scale. 


From the point P;, at the instant m — m; an 
observer is supposed to travel towards P,;, and at the 
instant m — mj,, another, also in the direction of 
P; from P;,,.'The two will arrive at P; at the instant m. 
The straight line g; corresponds to the observer 
moving from P;, in the positive direction of x, 
whereas the line g_ corresponds to an _ observer 
travelling in the opposite direction. Since the two 
observers, on arriving in P;, both experience the 
same current and voltage, the co-ordinates of the 
point of intersection of g,; and g_ are i and u, the 


current and voltage at P; at the instant m. 


w;(m) will now be taken to symbolize the state 
of the line at the point P; at the instant m. The 
symbol stands for both the numbers i(.x,;,m) and 
u(x;, m). According to the remarks above, the follow- 
ing principle applies: 

When -4(m — m) and yw 4,(m — m+) are 
known, \;(m) is determined by a simple algebraic 
operation, the procedure for which is as follows. 


@_: w— Zi=const. 


ball —To 1) 


ap;-1(t-7;) 


&4:uU+ Zt=const. 


BROWN BOVERI 


117277+1 


Fig. 1. — Explaining the Schnyder—Bergeron method 


If an observer were to move along the line in the positive direc- 


tion at the speed of propagation c, he would at all times experi- 


ence the relationship between the voltage u and current i as 
expressed by the line g;. Accordingly the line g- would apply 
for a similar observer moving in the opposite direction. 


Z = Line impedance 
From y;-,(m — m;) and ;,,(m — m;,,) we can 


determine the two straight lines g, and g_: 


For g; a,i+tbi,u+e,=0 


Forg. a_i+b_u+c_=0 


The coefficients a,, etc., are constants whose values 
can easily be calculated; they are governed by Z and 
the states y;,(m — m;) and ¥;,,(m — m4) (cf 
Fig. 1). Resolved for 7 and u these equations provide 
the new state y; (m), which can be written sym- 
bolically as 
vj (m) = S (g-, 8+) 

In which § is a symbol for calculating the point of 
intersection of the two straight lines. The two ends 
P, and P,, must be dealt with separately because here 
boundary conditions must be taken into account. 


In general form these are introduced as follows: 
q; (i, u,m) = 0 fa Oe 


in which q; is a function of the three variables 


1, u and m.+ 


4 In the largely automatic program described in section 4, 
the form A,(m)i + A,(m)u + A3(m) = 0 is assumed for q. 


P, P, | ne Sie SE Pa 
+ tt = ee 
Xo x Kim Kf Apes Fae 

q T Tj+r —— x Te, 
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Fig. 2. — The line examined 


P; (where j = 0...n) denotes the point on the line having the co- 
ordinate x;, and tj (where j = 1...n) the transit time of the 
wave between P;-, and P;. 
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Fig. 3.-— General arrangement of the line covered by the program 


a= Quadripoles consisting of simple basic networks 
L= Line or cable 
b =Symbol for boundary conditions 


These boundary conditions include such condi- 
tions as short circuit, line at no-load and connection 
of the line to an infinite bus at voltage w. 

To determine the current and voltage at Po, the 
following procedure is adopted. 

For an observer moving from P, towards P, at the 
instant m — m, a straight line g_ applies. He arrives 
at Py at the time m, where the relationship qp (2, u, m) 
applies. Hence the operation 


S (@o; 8-) 
yields the current and voltage at the time m. This 
time the S symbol is used in a somewhat extended 
sense. Analogue considerations apply to P,,. 


3.3 Expansion of the Theory to General Arrangements 


Fig. 3 shows the general diagram dealt with by 
the first program to be set up for problems of this 
kind. Line sections with generally different impe- 
dances are connected to one another by means of 
certain simple basic circuits. These circuits are 
referred to as quadripoles. In this general case too, 
the calculation of the conditions can be referred back 
to the S operation dealt with in the foregoing section. 

The suffix j is now used to number the various 
junctions, each of which is allocated two states 
Wi,;(m) and w_,;(m), as shown in Fig. 4. An 
ordinary point P; on the line is dealt with as though 
the line sections to left and right of P; were joined 


J 


by an “empty” quadripole with direct-axis_re- 


actance = 0. 


117280-1 
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Fig. 4. — Notation of a quadripole 


The functions w-, ; and w +, ; stand as symbols for the corre- 
sponding currents and voltages. 


4 


To calculate y,,; a modified straight line g* is 
used which, together with g, determines the state 
W4, ; (m) by means of the S operation. w_, ; (m) can 
finally be obtained by simple algebraic operations. 

If the quadripole consists of straightforward re- 
sistances, the coefficients of g* are easily determined 
from g, and the coefficients of the quadripole, for 
which g_ represents an unmodified straight line 
used for the S operation at an ordinary point on the 
line. In this case the calculation is strictly correct. 

If, on the other hand, the quadripole contains 
capacitances and inductances, the use of a modified 
straight line g* is only correct as an approximation. 
In this case the coefficients of this line, in addition 
to the data of the quadripole and the coefficients 
of g,, also contains the properties of state (currents, 
voltages) of the quadripole at the instant ¢ — Tp, if 
the state at time ¢ has to be calculated. Here Tp is the 
elementary transit time introduced in 3.2; thus all 
transit times t; are whole-number multiples of Tp. 
The smaller t), the better the accuracy of the ap- 
proximation. The most suitable value for tT, will in 
most cases be determined by the quadripole with 
capacitances and inductances; regarding the difficult 
question of the correct choice of Ty some hints are 


given in section 5, with reference to examples. 


3.4 The Block Diagram 


Let us suppose m,,,, 18 the largest of the numbers m; 


max 
defined in 3.2. The switching action which initiates 
the transient takes place at ¢ = 0. It is assumed that 
the states y,,; (m) and y_,; (m) are known for 


m= —m .. — 1. The state of the line corre- 


max * 
sponding to the longest transit time m,,,, T) must 
therefore be known too. This is fulfilled when the 
transient follows a steady state. This assumption 
regarding the past history of the line may appear 
rather odd, since it ought to be sufficient to stipulate 
that the initial condition at t = 0 be known. However 
it must be borne in mind that, in accordance with 
the assumptions, only the states at the point P; are 
known and not those of the whole line at the time 
t = 0; the states before that time provide the neces- 
sary additional information. 

If the now known states in the past are imagined to 
be arranged in a rectangular matrix, in rows with 


respect to time in the correct chronological order, 
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the uppermost, empty line before the start of the 
calculation being allocated to t = 0 and the bottom 
line to f= —T)™m,,,,. Beginning at one end of the 
line, the states at ¢ = 0 following the switching 
operation can now be calculated, and the upper- 
most line thus filled out. Having displaced each row 
downwards by one position, as a result of which 
the bottom row can be omitted as superfluous, the 
top row can be calculated again in the same manner, 
thereby yielding the states at the time ¢ = 1 . Tp. 
Repetition of this process finally yields the states at 
any desired instant t = mT ). Naturally in the com- 
puter the matrix corresponds to a definite storage 
sequence, for which additional storage cells have to 
be reserved for the properties of state of the quadri- 


pole with capacitances and inductances referred to 
leh sec 


4. Short Description of the Program 


The organization of the program is such that the 
cases occurring most frequently can be dealt with 
without additional programming being necessary, 
thus restricting the preparatory work to the produc- 
tion of a punched tape for the input data. Under 
certain circumstances it may be necessary to prepare 
special sub-routines for the boundary conditions. 

The topological structure of the arrangements 
which can be handled corresponds to Fig. 3, in other 
words only lines with no branches can be dealt with. 
Up to 24 line sections with different impedances 
and different transit times bearing definite ratios to 
one another can be handled. Between each line sec- 
tion a four-wire junction is visualized. Not more than 
nine different types may occur, though the same type 
may of course repeatedly occur, with different con- 
stants. Combined on a sub-routine tape are the 
calculation instructions for the modification of the 
straight lines g, as mentioned in 3.3, for nine typical 
quadripoles. Each type is allotted a reference number 
from 1 to 9. 

A punched tape is now produced, containing the 
following information: Transit times and impedances 
of the various lines, boundary and initial conditions, 
print-out characteristic. The quadripoles are de- 
fined by their data and computation instructions. 
The latter determine the sequence and types of the 


ake eS St 
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Fig. 5. — Representation of the examined line by a characteristic 


The above diagram represents the characteristic | 1 151116. 


lL 1.1 5a 
denotes an arrangement of the form illustrated 


various quadripoles; for example, 
in Fig. 5. 


5. Examples 


General Remarks regarding the Examples 


The examples quoted were primarily used for 
checking the program and the method. In this case 
it is not enough merely to check the program for 
logical correctness. Most examples contain quadri- 
poles for which approximation methods were used, 
as mentioned in 3.3. For such programs it is always 
essential to check through a few test cases, during 
which the necessary experience is gained, and con- 
fidence in the method. 


Line Terminated by an Inductance (Fig. 6) 


The line terminated by an inductance LZ and dead 
att < 0, hasaconstant voltageu = 1 applied att = 0. 
The current 7 through the inductance is to be cal- 
culated. This problem can be solved relatively easily 
with the method using eigen functions; the numerical 
solution is nevertheless so complicated that it is 
worth programming. The current through the induc- 
tance was deliberately calculated instead of the 
voltage, because the voltage, being a discontinuous 
function, as described in section 2, cannot be satis- 


factorily expanded with eigen functions. 


The calculation was based on the following 


numerical values: 


i L 


/ 


; 
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Fig. 6. — Check calculation to verify the program 


At ¢ = 0 the constant voltage U is applied to the line. The 


current in the inductance L is calculated. 
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Fig. 7. — Result of the check calculation as in Fig. 6 


Current i in function of time ¢. The time unit is t,, the transit 

time of the line. Calculation by the Schnyder-Bergeron method 

and the development of the solution using eigen functions gave 
approximately the same result. 


Line length lkm  (t,= 3-33 x 10° s) 


Impedance 500 ohms 
Inductance L 0.0167 H 
Applied voltage U 1V 


Fig. 7 shows the curve of the current in the reactor 
as a function of time, the transit time t, being chosen 
as time-scale. Up to 31, the calculation can be 


checked direct, since on expiry of one transit period 
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Fig. 9. — An example used to check the program 


The line is terminated by a resonant network tuned to the power 

frequency. Parallel to the latter is a resistance, whose value is 

equal to that of the impedance. The line is connected to the a.c. 

source U cos wt. When the transient has decayed, the line trans- 
mits its natural power. 


Length of line 750 km Impedance Z = 600 ohms 


T,, during which the current remains zero, it rises 
according to the well-known theory, in accordance 
with the law [8] 


vA 
i= (inet | 


This law remains valied until the wave, reflected 
at t = 7, at the voltage source, returns for the first 
time with the opposite sign to the inductance L. 


Us ; 
Since the time constant Z = 3-33 & 10 s is 


large compared with t,, the current rises almost 
linearly. The rest of the curve can also be explained 
fairly easily. 

Obviously the ratio of ty to the time constant L/Z 
governs the accuracy of the approximation method 
described in 3.3. In the present case it was made very 
small, i.e. 0-002. Increasing it to 0-01 produced slight 


differences which are just visible in the trace. 


Fig. 8. — Approximate method of 
taking line losses into account 0 


Showing the current at the infeed 


point when a 600 km line with losses 
is connected up. Due to the damping 
caused by the ohmic line losses, 


the discontinuity becomes gradually 
smaller. 

Z=650 ohms R=0-104 ohm/km 30° 1 

u=1coswt w= 314 s-} 
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Fig. 10. — Curves showing the 
currents i, and iy in the circuit in — 
Fig. 9, as a function of the time t — 
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i 

ae 
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In Fig. 10a it can be clearly seen — 
that the current, which is severely 
distorted to begin with, gradually : 

‘ 


=| reverts to the pure sine wave. (Small 
capacitance C = 0-203 x 10-6 F.) 

In Fig. 10b hardly any damping” 

3 ats 4 —L can be detected. (Large capacitance 


0 ] 2 3 5 6 7 S 9 10. -10°s 
BROWN BOVERI a. aes 117302+1 C = 0:203.x 104.8} F 


Line with Losses 


Line losses can be taken into account by dividing on a 600 km long line proved that division into only 5 
the line and connecting a resistance as a quadripole _ sections yielded almost the same results as a finer 
between each two sections. Calculations carried out division into 15 sections. 


= : 
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Fig. 11.— Investigation of transients 


following restriking in the breaker S ; 
tC; 


Before the restrike the two line sec- 
tions P, P, and P,P; are equally and 
oppositely charged. For notation, see 


u= cos wl 


text. 
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Fig. 8 shows the current produced when a 600 km 
length of line at no-load is connected to a 50-c/s 
supply at the instant when the voltage is at a maxi- 
mum. It will be observed that discontinuity de- 
creases with time, the curve approaching gradually 


to a sine curve. 


Application of an A.C. Voltage to a Line 
Terminated with a Resonant Network 


The arrangement is illustrated in Fig. 9. The 
resonance frequency of the termination is equal to 
the applied frequency of 50 c/s and the resistance is 
equal to the impedance. This theoretical example 
was employed for check purposes because the oscil- 
latory state is very simple. As can be easily seen, it 
corresponds to the transmission of the natural power. 
The transient is damped by the resistance in the 
terminating network, in this case the impedance. 
During the calculation it was confirmed that the 
transients only decay in a short time, due to the 
damping effect, when the capacitance is sufficiently 
small. The damping effect of large capacitors is very 
slight because they absorb the incoming wave and 
very little energy is dissipated in the parallel resis- 
tance. Fig. 10a shows the currents at the ends of the 
line when the capacitance is relatively small. It can 
be clearly seen that the very discontinuous curve at 
the beginning is damped within a few cycles to such 
an extent that only the steady-state curve remains. 

The step length t, was made 10~ s. In this case 
the time constant, formed by the impedance Z and 
capacitance C, governed the estimation of tT). (The 
time constant of the inductance is much larger in 
this case.) Since C = 0-203 x 10-6 F, and Z = 600 
ohms, T= 1-218 x 10-*s. Therefore t) and T; are 
both of the same order of magnitude. 

By way of comparison, Fig. 10b shows the current 
i, at the line termination when the capacitance C is 
made 100 times larger. Here, for the reasons given 


117287*1 


above, the effect of the damping is almost imper- 
ceptible. 


The Effect of Restriking in a Circuit -Breaker 


According to Fig. 11] it is assumed that at t = 0 
the breaker S§ is closed, the line sections on either 
side of the breaker being at equal and opposite 
potentials before closure. 


Data: Impedance Z = 500 ohms 
Line length /,; = 1 km 
l= Lkm 


Capacitance C, = 800 pF = 8 x 10°7°°F 
Go 142 oF = 142% 104 F 
Inductance L, = 0-36 H 


Voltages before closure of the breaker: ° 
Section PoP, +1V 
PEP eV 


Fig. 12a to e show the curves of the current 
and voltage at the various points in the diagram. 
Again it can be established that, agreeing with the 
theory, the incidence of a wave at the capacitors 
results in an exponential shape of the current and 
voltage curves. Moreover, at P, the inductance is 
so high that, to begin with, it cannot affect this shape. 
The curves shown dotted in Fig. 12a—d correspond 
to the case in which the line is open at P, (small 
capacitance and high inductance) and _ short-cir- 
cuited at P, (large capacitance). This case, which 
provides a very clear picture, makes it much easier 
to understand the complicated phenomena. For this 
example the choice of the basic transit time had to 
take into account the time constant T, = C,Z, 
which is very small. The ratio t): T, = 0-139 gave 


satisfactory results. 


® Since the circuit is linear, the voltage level of the initial 


conditions can be selected arbitrarily. 
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Fig. 12. — Transients caused by restriking 


: uy = Voltage at P, 
: Ug = Voltage at Py 
: t,= Current at P, 
iz = Current at P; 


como OT DP 


: ic, = Current in C, 


The time is shown as multiples of t, = T., the transit times of 
the two lines. The curves shown dotted indicate what happens 
when the line is open at P, and short-circuited at P3. 
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6. Extension of the Program 


Since the present article was first drafted, the 
program described in section 4 has been consider- 
ably extended. The conditions stipulated regarding 
the topology of the network, as shown in Fig. 3, were 
discarded. Now it is possible to investigate systems 
with any number of branches, as illustrated by a 
simple example in Fig. 13. 

The original intention was to utilize the program 
for the investigation of overvoltage phenomena in 
connection with line protection by lightning ar- 
resters. On this account it was necessary to make a 
further sub-routine capable of reproducing the 
properties of an arrester with reasonable accuracy. 
It is well known that an arrester acts as an insulator 
until the moment when the breakdown voltage is 
reached; subsequently the voltage is determined by 
the current flowing through the arrester. The dis- 
charge characteristic is represented in the program 
by a succession of straight lines, as shown in Fig. 15b, 
for instance. 

The principle of the extended program will now 
be explained by reference to an example. The layout 
is illustrated in Fig. 14. At point P, there is a 220-kV 
transformer, represented by an input capacitance 
C and impedance R. At a distance of 15 m from the 
transformer, at point P,, there is a branch leading 


Fig. 14. — Diagram for the investi- 
gation of overvoltage phenomena in 
a 220-kV installation 


BROWN BOVERI 


Point P;: Connection to a transfor- 
mer 
and R 


represented by C 


Point P,: Branch 
arrester with non-linear 
resistance R (i) 


to a_ lightning 


Point P,;: Connection to transmis- a: Z,=Z,=Z3=Z, = 400 ohms 
Cy = Cg= Cg= C= 300 m/s 
%=7,=0°05 ps t3=2 us 


sion lines represented by 
the resistance R, 


a 
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Point Py: Point treated as the start 
of the line in the event of 
a lightning stroke 
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Fig. 13. — Example showing the topological structure of the 
system which can be examined with the aid of the extended 
program 


i= Line or cable 


a= Electrical circuits or boundary conditions 


to the lightning arrester 15 m away. 90 m before 
this branch, at P,, are the connections to very long 
transmission lines. These can be represented by a 
resistance A, corresponding to their impedance. 
The cases investigated were R, = c (no branches) 
and R, = 200 ohms (two parallel lines). At a point 
600 m from P, lightning is assumed to strike the 
line at Py, the amplitude being limited by the insula- 
tion level of the line. This was assumed to have the 
very high value of 1500 kV and the rate of rise of 
the incident wave-front 1000 kV/us, a value com- 
monly used for calculations when unfavourable con- 
ditions are anticipated. The shape of the incident 
wave is illustrated in Fig. 15a, while 15b shows the 


characteristic of the arrester. It was assumed that 
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T, =0°3 us; R=7 kohm 
C= 1000 pF 


R,=c 


Two examples a and b were calcu- 
lated, using the following numerical 
values: 


For the fundamental shape of the 

voltage u(t) see Fig. 15a; for the as- 

sumed arrester characteristic u(i) see 
Fig. 15b. 


b: As for a, but with R,; = 200 ohms 
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Fig. 15a. — Assumed shape of the voltage u(t) 


the latter has a protective rating of 680 kV, con- 
forming to the Swiss co-ordination regulations for 
220-kV installations with effective earthing. Thus 
at the rated discharge capacity of the arrester, 10 kA, 
the residual voltage should just amount to 680 kV; 
the actual breakdown voltage is 660 kV. The quan- 
tities calculated were the voltages at the transformer 
(uy) and arrester (w,) and the current through the 
arrester (7,). The results are plotted in Fig. 16, I and 
II. It will be observed in I that the voltage at the 
arrester never exceeds the protective rating. Since 
the transformer is 30 m away from the arrester, the 
voltage uy rises briefly to 910 kV. According to a 
commonly used rule-of-thumb, for a distance of 
30 m and a rate of rise of 1000 kV/us, a brief voltage 
30 
300 


pected. The calculation, however, proves that this 


l 
rise of 2 . 1000 = 2 1000 = 200 kV may be ex- 


LOOOr1L0 


400 


200 


kA 10 
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Fig. 15b. — Assumed characteristic u(t) of the arrester 


rise may actually be somewhat higher. The illustrated 


eS. Sa a 


7 


voltage curve at the transformer corresponds to a_ 


full wave with an amplitude of about 800 kV. From 
Fig. 16 II it is possible to detect the influence of the 
parallel lines connected at P,. The rate of rise of the 
incident wave is reduced by half. The result of this 
is that the maximum voltage experienced by the 
transformer is only about 785 kV. One might per- 
haps tend to expect a still greater influence from 
the parallel lines. The reason why this is not so is 


that the amplitude of the incident wave was as-— 


sumed to be roughly twice as high as the protective 
rating of the arrester; hence two parallel lines are 
unable to reduce the voltages to such an extent that 


the arrester does not have to function. 


kV|kA 


se, se 


pols 


BROWN BOVERI 


1000r10 


3 


7293-1 


kV|kA 


curves at the transformer (ur), 
at point P, (uy) and at the arrester 
(u,), as well as the current i, 

through the arrester 


= 
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3 I: for example a in Fig. 14 
II: for example b in Fig. 14 
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Fig. 16. — Calculated voltage/time 
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7. Conclusions 


Investigations carried out have confirmed that 
for discontinuous phenomena the method employed 
here is more suitable for the particular problems 
dealt with than the expansion of eigen functions. If 
the circuits contain capacitance or inductance, the 
Schnyder—Bergeron method can only be employed 
in practice on a digital computer because, in view of 
the accuracy required, the calculation has to be 
split into a very large number of steps. For this 
reason the time taken to complete some of the calcula- 
tions listed in the article as examples was up to one 
hour. But it has been shown that a method origin- 
ally intended as a graphical solution can be suc- 
cessfully adapted to render it applicable to a digital 
computer. 

W. FREY 


P. ALTHAMMER 


(KME) 
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INVESTIGATING THE TRANSIENT STABILITY OF SYNCHRONOUS 
MACHINES WITH THE AID OF A DIGITAL COMPUTER 


The article describes a program for calculating the transient 
stability of a synchronous machine connected through a 
general quadripole to an infinite bus. To illustrate applica- 
tions, examples are quoted of common methods of calculating 
the stability and of the stability with automatic rapid 


reclosure. 


1. Treatment of Stability Problems 


HE CALCULATION of the behaviour of syn- 
chronous machines in the face of disturbances 
can only be performed in a reasonable length of time 
with the aid of electronic computers. Only isolated 
cases, involving say two or three machines, can be 
tackled rationally without such modern aids. But 
even in these simple cases quite extensive simplifica- 
tions have to be made. 

To determine the transient stability of a system 
with several machines the a.c. network analyser 
described in the article beginning on page 313 has 
given good results for some time, but lately the digital 
computer has begun to offer serious competition. 
The investigations on the analyser are also subject 
to quite considerable simplifications. For example, 
owing to these simplifications, it is impossible to 
determine the effect of the control system, the damp- 
ing winding and even of the synchronous reactances 
on the transient stability. 

The accuracy attained with the network analyser 
is generally adequate when the stability of a whole 
network has to be investigated. In most cases no 
higher accuracy could be expected from refined 
methods, because when whole networks are being 
examined some of the essential data are not known 
to any great accuracy. 

However, in many cases an improved method for 
calculation of the stability is desirable, particularly 


681.14: 621.313.32.016.352 


when the influence of certain parameters on the 
stability has to be determined. In such problems it 
is often only a single machine or power station which _ 
is connected to an infinite bus by means of a trans- 
mission line. In such cases an ideal means is the 
analogue computer [1]. In certain cases, however, 
setting up the problem and obtaining the measure-— 
ments on the analogue computer prove rather 
tedious and long-winded. A case in point is the prob-_ 
lem of automatic rapid reclosure, where three | 
different states have to be considered: normal 
operation, short circuit, interruption. If the effect of 
different parameters has also to be studied—in-— 
volving the performance of a large number of calcu- | 
lations each time—the inautomatic operation of the 
analogue computer proves quite a serious handicap. — 
For these reasons a program for investigation of | 
the transient behaviour of a single synchronous j 
machine has been prepared. In contrast to the - 
analogue computer, once a program has _ been 
checked, it does not occupy any further preparation — 
time at the digital computer, from which arises the f 


important advantage that relevant enquiries can ~ 
be answered in quite a short time. ) 
The mathematical fundamentals of the method 


and program will now be described, reference being 
made to the literature for certain details. 


2. Stipulations for the Calculation 

The arrangement handled by the program is 
illustrated in Fig. 1: A synchronous machine S is 
connected with an infinite bus N through a pas- 
sive network represented by V“. The voltage regu- 
lator and the exciter are regarded as being combined 
in the excitation circuit E. The superscript 7 in the 


¢ 
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ig. 1. — Schematic diagram to illustrate the set-up dealt with 
by the program 


S§ = Synchronous machine 
E=Excitation circuit, comprising exciter and regulator 
V@ =Quadripole representing the connection between the 
machine and the network with an infinite bus 
N = Infinite bus 
u, = Terminal voltage 
us) = Voltage of the infinite bus during the 7-th state 


i=Superscript denoting the system state 


The quadripole V() may, for instance, represent a transformer 
with one or more lines. A fault in the system, e.g. interruption 
of a line, can be represented by an appropriate change in the 
quadripole. Loads can likewise be represented by the quadripole. 


case of V and u®) indicates that, during the tran- 
‘sient phenomenon, both the quadripole and _ the 
system voltage may vary. Henceforth the various 
states are referred to as system states, the state im- 
mediately preceding the disturbance being denoted 
by the superscript 0. When auto-reclosure is suc- 
cessfully executed, three different system states are 
jinvolved, namely V“) before the disturbance, V“) 
during the short circuit, V® during the process of 
jinterruption; following reclosure V“) is once more 
festablished. In the event of a persistent fault the 


}sequence is V — VO — V® — VO — Ve, 


2.1 The Synchronous Machine 


The mathematical representation of the syn- 
chronous machine is in accordance with the well- 
known Park’s theory, the equations of which can, 
however, be considerably simplified for this pur- 


pose. Here reference is expressly made to [3], in 
which the equations and notation are given. As will 
| be known, this theory also takes into account the 
: effect of the damping windings in the two axes, each 
} of these windings being represented by a circuit of its 
wn. The saturation of the machine, on the other 
and, is not allowed for, as is also the case in the 


present program. However, the capacity of the 
digital computer employed is large enough to per- 


mit extension of the program to include saturation. 


To characterize the synchronous machine the 


following constants were fed into the computer: 


Xd Xq = synchronous reactances 


x’), x"; = transient and subtransient reactances in 


the direct axis 
x", =quadrature-axis subtransient reactance 
x, = leakage reactance of the stator winding 


and 


short-circuit time constants 


T’»,, T’,=direct-axis transient subtransient 
T”, = quadrature-axis subtransient short-circuit 
time constant 


T,,=starting time constant of the machine set 


The above are known as the operational con- 
stants: apart from x,, x, and T,, these constants are 
derived quantities which do not appear in the funda- 
mental equations of the machine. For the sake of 
brevity the constants appearing in the fundamental 
equations will be called basic constants. These are 
the various leakage coefficients, x,, x,, and the time 
constants of the different windings, when the other 
circuits coupled to the particular winding are open- 
circuited. 

As was demonstrated elsewhere [2], the basic 
constants can be determined back from the opera- 
tional constants, though for this an insignificant 
supplementary hypothesis has to be made. The 
resistance in the rotor of the machine is also neglected. 
If necessary, however, it can be taken into account 
in the quadripoles V“. The program begins with 


the computation of the basic constants. 


2.2 The Excitation Circuit 


The parts of the program concerning the excita- 
tion are incorporated in the main program as sub- 
routines, enabling a variety of excitation systems to 
be tried out, if required. In the examples quoted in 
the present article a regulator with proportional- 
integral action (PI), whose measuring system has a 


1 This is not fulfilled with the operational time constants 
Te etc. 
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time constant, and an unsaturated exciter were 
assumed. The deflection of the regulator is limited 


in both directions. 


2.3 The Quadripole V 


The quadripole V can be characterized by the 
four components of its iterative matrix (see page 359) 
which are complex in the general case; their depend- 
ence on frequency is not taken into consideration. 
Accordingly, for each system state, four complex 
numbers have to be stored. In order to give the 
program the necessary flexibility, the voltage of the 
infinite bus is also stored along with the components; 


this voltage may also be subject to certain fluctuations. 


2.4 The Mathematical Method 


The complete synchronous machine with damping 
windings in both axes leads to a differential equation 
of the 5th order. An assumption made regarding 
the excitation circuit was that it could be described 
by a system not exceeding the 10th order. The over- 
all system of not higher than the 15th order is inte- 
grated step by step using Runge-Kutta formulae of 
the 4th order. 


3. The System of Equations 


3.1 The System for the Synchronous Machine 


To formulate the Park system of equations, prior 
knowledge of which is taken for granted, the follow- 


ing notation is introduced: 


Winding or axis: Subscript? 

d axis of the stator | D; 
Field winding uae ma 
Damping winding in the direct axis 3 aa 


q axis of the stator 4 
Quadrature 


Damping winding in the : 
| axis 


quadrature axis 5 


1, where j= 1.90), Currents 


jJ=1...5 Voltages (u, = u; = 0) 


Uu; ” 


Wis * j=1...5 Flux linkages 


* To economize in the notation it is preferable to employ 
numbers instead of letters as subscripts. 


T;, where j=2,3,5 ‘Time constants in the rotor 

windings 

T, Starting time constant of the 
machine 

m, Driving torque : 

w, Rated angular frequency 

6 Angle between the induced | 

voltage and the fixed system | 
voltage 

uy System voltage’ 

u, ‘Terminal voltage 


ly, where j,k=1...5 Induction coefficients 


The induction coefficients are zero when j and k refi 
to two different axes, i.e. when j < 3 and k > 4, and 
vice versa. Apart from the time constants, all values. 
are expressed in the per unit system (for details see’ 
[3] ). The coefficients J, are either basic constants 
or simple combinations of these. 


Employing this notation the system of equations: 


&. 


for the machine can be written as follows: 


= | dy, Ya ( do | 

Uy = Somr a, On or =) al 

a l dy, vy awd | 

a (oo) Sar On (, a re F 

dy; : 

w= i+ Te, = 9,3, 5 (2) 

5 

v=) lea, ey lee (3) 
k=1 

, : Toe aaG 

My = Yay — Vila + Per. (4)) 


Remarks regarding the system (1) to (4) 


— Incontrast to reference [3] in the bibliography 
the time ¢ is not introduced in the per unit system 
but is expressed in seconds. For this reason m, appe 


sometimes in the denominator. 
: . dd 
— In [3] the equation (1) contains “aH 
d2 
ee instead of 
&® was defined as the angle between the fiel 


instead of «,, + eu and in (4) 


® Vectors of current or voltage are denoted by an 
above the symbol. 
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winding and the winding axis of one of the three 
phases of the stator winding, wherein, without im- 
posing any restriction on the generality, a two-pole 
machine is assumed. Therefore 


IU 


=> w,t + 6 oy ’ 
from which 
di do d?p d*6 
Seemed ed an 
— In (1) the terms a and ~ describe 


the transformatory effect on the stator windings of 
changes in current in the rotor. In stability investi- 
gations it is permissible to suppress these terms. 
Likewise the derivatives of 6 with respect to time may 
be neglected compared with «,. 


Thus (1) simplifies down to 


ahi het | (1a) 
Sat 


“= 


According to the vector diagram (Fig. 2) 


u, = uy sind,, | (5) 


U, = uy, cosd 


m 


in which 0,, is the phase angle between the induced 


e.m.f. ¢ and the terminal voltage. If the machine 
is connected direct to an infinite bus, i.e. when uy 
is constant, the equations (la) and (2) to (5) are 
sufficient for calculating the phenomena in the 
synchronous machine. In this case 6,,= 6. This 
assumes that the voltage u, applied to the field 
winding and the driving torque m, are given. If 
tt, is not constant, (5) must be replaced by a system 
containing the properties of the quadripole V“. 


3.2 The Equations of the Quadripole V 


Using the notation in Fig. 3 the following equations 
apply for the quadripole: 
Sas Ay uy + Ap iy 
ly = ae un a ae iy 
The bars denote that these quantities are of a com- 
plex nature. 
eee An 
Ay» 


Bs | 


22 


jixa 
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Fig. 2. — Vector diagram of the synchronous machine 


In order to simplify the method of writing the equations, the 
voltages in the two axes are denoted by wu, and uw, instead of 
the usual uy and ug. 


it follows that 

“lat eAles, — i 

ue = Dy uy + Dy t 
having regard to the fact that the determinant of 
the iterative matrix is unity. 

This vectorial equation independent of the sys- 
tem of coordinates is resolved into two scalar 
equations. In the system of coordinates whose axes 
coincide with those of the machine, with the above 


notation, it follows that: 


Dy =o + Jé, 


Dy =1 + jx 
u, = ry — xt, + (@ sind — &cos 6) uy | 
Uy = xl, + rig + (0 cos 6 + € sin 0) uy 


In the general case the equations (6) supersede 
those in (5). 


3.3 Integration of the System 


In the system (la), (2), (3), (4), (6) only the 
equations (2) and (4) are differential equations. Let 


us substitute 


117240 | 


BROWN BOVERI 


Fig. 3. Notation of the quadripole 
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in which v is the speed deviation during transient 
phenomena referred to the synchronous speed. ‘Thus 
the differential equations of the system can be ex- 


pressed in the form needed for integration 
dy; 
dt 


du 
dt 


sa 
di 


wy (a; 4) 59 = 2,'35.9 


= Ti} (my — Yat + Yi 4s) (7) 


— WO, U 


The main variables 


Of the various variables of the complete system, 
five, i.e. the three rotor fluxes, the angle 6 and the 
speed deviation v, are specially distinguished, in that 
only derivatives of these variables occur. They will 
therefore be referred to henceforth as the main 
variables. 

For the computation procedure it is now most 
important that, when the state of the system changes, 
these ma‘n variables do not change abruptly like the 
others. This can be seen from the above equations, 
which retain their validity when the state of the 


system suddenly changes. Since their right-hand 


117241-1 


Fig. 4. — Flow diagram for calculation of the right-hand sides 
of the equations (7) 
The quantities in the blocks are calculated in the indicated 


part of the program. The numerals refer to the relevant equa- 
tions in the text. 


From the organization aspect, this program is a sub-routine 
of the Runge-Kutta sub-routine. 


I: from the Runge-Kutta sub-routine 
II: return to Runge-Kutta sub-routine 


sides remain finite, this excludes any abrupt change 
in the variables. For the quantities y; this dis- 
covery complies with the principle of constant — 
flux linkage. Since the derivatives of y, and y4 are— 
not taken into account in equation (1), the stator 
fluxes can, on the other hand, change abruptly. 

The constancy of the main variables when the 
state of the system changes, i.e. corresponding to a 
change in the system of differential equations, 
renders it unnecessary to re-calculate the initial 
conditions after the change. 

The main variables determine the state of the 
machine completely at all times. Only the main 
variables are stored in the program for the Runge- 
Kutta method. i 


Starting with these variables the right-hand sides _ 
of the system (7) are calculated at each step, ac- 
cording to the following procedure: 


currents 7;: 
Ti S Aig Ve (8) 
k 


} 
, 
The system (3) can be resolved in terms of the 
Since the matrix of the coefficients /, contains a 
large number of noughts, simple expressions can be . 
used for /;, in the inverse matrix. : 
In (6) the substitution corresponding to (la) is | 
made for the left-hand sides, while 7, and7, are — 
replaced on the right-hand side by the fluxes ,, with . 
i 


the aid of equation (8). 


We obtain a system of the form i 


Ay, Wy + Aq Ya = Q | 
Ag, Yy + Ogg Ya = Q | 
with 


Q1 = Me Wo + Ys Ys + Us Ys + 
+ euysind — €uycosd 


(10) 
Qs = M2 V2 + 443 Y3 + a5 Y's — 
— ouycosd — Euysind 
The a,x bear a simple relationship to A, as well as 
to r, x, @ and é. 


Resolving the system (9) for y, and y4 gives the 
two stator fluxes as a function of the main variables. — 
Inserting these values in (8) yields the currents 7;, so” 
that now all variables have been determined as a 
function of the main variables. With these the right 


hand sides of the equations (7) can now be calculated 
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Fig. 4 shows the flow diagram of the rather com- 
plicated computation of the right-hand side of (7), 
denoted in the diagram by /;. The necessary con- 
stants, such as aj, etc., were calculated before com- 
putation started. 

The general Runge-Kutta sub-routine used _ is 
based on the following assumptions: 


The system to be integrated is of the form 
Ss Pe we ees 2 


The following information must be inserted: 


Number of equations n 
Length of step h 


Initial addresses of the two series 


of variables x; and Si 


Address of a sub-routine calculating /; 


from x; each time 


When run through once the Runge-Kutta sub- 
routine, starting from the values x; at time ¢, calcu- 


lates the x; at time ¢ +h, the sub-routine being 


run through four times to calculate f;. 


3.4 System of Equations of the Excitation Circuit 


The terminal voltage u, needed in the event of 


voltage regulation works out in its absolute value to 


(see equation | a) 


u, = Vu + uj = Vw +y? 


The control system with PI action utilized in subse- 
quent examples gives three additional differential 
equations, one each for the time-lags of the measuring 
system, the exciter and the regulator with PI action. 
The equations are simple and do not need to be 
listed here. Programming the limits merely needs 
a little thought. 


3.5 Initial Conditions 


The initial values of the main variables and the 
additional variables are computed by the program 


before integration starts. 
‘ 
For this, with the other data the following values 


be given, referring to the state before the dis- 
nee (superscript 0): 
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Fig. 5. — Flow diagram for calculation of the transient stability 


Various details, such as printing instructions, have been omitted; 
likewise the links depending on whether a damping winding is 
provided or not. 


po = active power 
uy = voltage of the infinite bus 


uw) = terminal voltage of the machine 
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From the machine data and the quadripole co- 
efficients A) of the quadripole V we can calcu- 
late all the variables needed. Here it may be pointed 
out that, mathematically, the calculation yields two 
solutions; of these we have to choose the one corre- 
sponding to stable steady-state operation, because 


at least one of them is unstable in the steady state. 


4. The Flow Diagram 


Omitting details, Fig. 5 shows the flow diagram 


for the entire program. 


4.1 General Remarks 


The various states of the system are characterized 
by the subscript 7, where 0 <7 < 10, to which a 
counting register is allotted. 7 = 0 corresponds to 
the state of the system before a fault or disturbance. 
i- 


1 


is the time at which the change from state 7 to 
state 7 + 1 takes place, and is included with iO in 
the input data; ¢, is the time up to which integi ation 


is to take place. 


4.2 Remarks on the Individual Blocks 


Block 1: ‘The basic constants are calculated from 


the operational constants (see section 2.1). 


0 05% Se ow kee ee ee Wire Ti Be 
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117243+1 
Fig. 6. — Short-circuit current during a three-phase short circuit 
when the synchronous machine is running at rated voltage 
The curve is the envelope of the symmetrical component. 


Machine data: xq = 1-6 per unit 


xd=0°3 p.u 
x= 02 p-u 
Ti=1-0s 
Ta=0-02 s 


Before the short circuit the machine was running idle at rated 
voltage, 


Block 2: The initial states are calculated from 
system state 0. 

Block 4: The counting register is set to 1, thereby 
inserting the first state after the distur- 
bance in the calculations. 

Block 6: The Runge-Kutta sub-routine | 

Block 6a: Sub-routine according to Fig. 4. . 


The blocks 2 and 6a also contain sub-routines for the — 


excitation circuit. : 


5. Examples 


5.1 Determining the Symmetrical Component 
of the Short-Circuit Current in the Event of a 
Three-Phase Short Circuit ; 

, 

For V the constants of a quadripole are inserted, ; 


given by a direct-axis reactance of x, = c (in 
practice it is 10"). V™ consists of a quadripole “a 
the quadrature-axis reactance xg = 0 (in practice 
10-!°). The sub-routine for the excitation circuit j 


corresponds to a constant voltage applied to the 


field of the synchronous machine, made equal to. 


1 per unit. With these assumptions the computer 


calculates the envelope of the symmetrical short- 
circuit current. The asymmetrical component is not © 
taken into account because, as explained in 3.1, the 
transformer effect is neglected. ; 

Fig. 6 shows the well-known curve of the short-— 
circuit current. This curve can be easily checked 
from the operational coefficients. | 

Calculation on the digital computer gave the 
following relative errors compared with calculation ~ 


according to the formulae: 


Time Relative error 
s 

0-04 0-02049 

1-00 ae Ns: 

4-00 < 10° 


Whereas for longer times the error is negligibly small, - 
at 0-04 s after the short circuit there was an error 
of 2%. The calculation was carried out with a stay 
length of 0-02 s, i.e. equal to the subtransient time) 
constant. The error determined for 0-04 s, as was’ 
also confirmed by later tests, is connected with 
rapid change in the current in the damping windi 


e obvious rule may be made that the length of 
e step should not be shorter than the subtransient 
ime constant. A remarkable fact is that an initial 
rror is not propagated over long periods. 


5.2 Examples showing the Procedure for 
Stability Calculations 


As mentioned at the start, simplifications are 
sually made when calculating the transient stability. 
Jonsidering a concrete example an investigation 
as undertaken to check the effect of this on the 
urve of the rotor angle with respect to time. 

It is assumed that a machine is connected to an 
finite bus by a two-circuit line 200 km long. A 
, ree-phase short circuit is assumed to occur at the 
iddle of one of the two lines. After 0-15 s the 
ffected line is locked out. The curves of the rotor 
ngle during this disturbance were calculated (Fig. 7). 
or the four curves the following conditions apply: 


: The system of equations must be complete. 


: The system of equations must be complete, but 
T", = T", = 0, thereby eliminating the effect of 
the damping winding. 

: The calculation is performed with constant flux 
linkage, without the influence of the damping 
winding. To obtain this variant we merely have 
Suaupstitute 7',— oc, T*,= 7", = 0. 
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4: Calculation conforms to the simplification com- 
monly used in practice, wherein only the transient 
reactance is effective in the two axes of the ma- 
chine. The variant was calculated with the 


following stipulations: 


,= T", =0, 1’, = ©, x, = *',,%; may have 


any value but it must be > x’). 


With this simplification the rotor angle is defined 
as the phase angle between the “air-gap voltage”’ 
following x’; and the fixed system voltage. During 
hunting it differs from the true rotor angle by 
a constant amount. This difference was not taken 


into account in Fig. 7. 


Data of the machine 


Machine rating = natural power capacity of the 


two-circuit line 


Synchronous reactances xy=1:0, x,=0-65 per 


unit 

Transient reactance Soy adr Dalle 

Subtransient reactances x" ,=0-2, x", =0-2 p.u. 
Leakage reactance of stator x,=0-17 p.u. 

Transient time constant ier les 

Subtransient time con- 

stants I", =0-02 s, T’,=0-02 s 


teen's 


Starting time constant 


ig. 7.— Rotor angle as a function 


of time 30 


-phase short circuit on one of 20 


0 parallel lines, cleared after 10 
15 s by disconnecting the affected 
line. 0 


Length of line 200 km. 


BROWN BOVERI 


tion assuming certain con- 
ditions. direct axis 
= Complete system of equations 


= Neglecting the damping winding 


3—=Constant flux linkage in the 


4 = Calculation with xin the direct 
and quadrature axis 


3,5 ee 


1172441 


x = Points calculated with slide- 
rule by the step-by-step 
method 
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Excitation circuit 


Transfer function = 


— Kae 


KS 10-7, = Ose Ue es 


] l ] ) 
ate: er ph ee 
PeaPpTy TEATS ( WET: 


Maximum exciter voltage + 5 p.u. 


— 1 p.u. 


Minimum - 9 


From curves | and 2 it will be observed that the 
damping winding only begins to take effect after 
one full oscillation, discernible in the increased 
damping of the oscillations. 

Curve 3, which for short times coincides with | 
and 2, shows that, to begin with, the machine 
exactly conforms to the principle of constant flux 
linkage. Of course, in this case the oscillations are 
undamped. Whereas in | and 2 the voltage regulator 
in time corrects the terminal voltage back to the 
rated value, this is not the case in 3 because in this 
calculation the effect of the regulator was ignored. 
This explains why the mean value of the rotor angle 
is smaller in 1 and 2 than it is in 3. 

Finally curve 4 shows the curve of the rotor angle, 
as obtained with the greatest simplification. It is 
surprising to note that the deflection of 6 is larger 
than in curve 3. It can, however, be shown that 
immediately following the short circuit the power 
output, when calculated according to condition 3, 
is about 27% higher than according to 4. 

The points (x) plotted in curve 4 were obtained 
by calculation with a slide-rule, using the normal, 
simple step-by-step method (step length 0-05 s). It is 
evident that this method yields surprisingly good 
results within the scope of the conditions stipulated. 
For calculations on the network analyser the same 


procedure is adopted. 


5.3 The Transient Stability Limit with 
Rapid Auto-reclosure 


5.3.1 Effect of the synchronous reactance on the transient 
stability in a special case 


The same conditions apply as in 5.2, i.e. the 
synchronous machine is connected to an infinite 
bus by a two-circuit line, the length of which is 
varied between 100 and 300 km, and a short circuit 


is assumed to occur at the middle of one line. But, 


—— 


in contrast to 5.2, it is assumed that after the line ha: 
been interrupted, the short circuit is cleared and th 
affected line can be reclosed after a dead time 0 


0-3 s. By using a program in which the origi 
program appears as a sub-routine it would i 
possible to work out the maximum power which ca 
be stably transmitted under the particular conditio 
Thereby the length of the line and the synchrono 
reactances x, and x, = 0°65 x,, were varied. Th 


remaining constants are the same as in 5.2. 


It must be pointed out that it does not entirel 
conform to the technical circumstances to vary thi 
synchronous reactance while keeping the transien 
reactance x’, fixed. However, in these calculation 
the main point was to determine the effect of th’ 
synchronous reactances on the transient stability 
For the same reason, powers were considered which 
exceeded the rated output of the machine. Fron 
Fig. 8 it can be seen that the transmissible powe 


P\ drops by about 15°% when the synchronov 
reactances are increased by a factor of 3-3, fror 
0-75 to 2:5 per unit. Although this dependence 
; 


upon since sometimes the view is expressed thay 


is not very marked, this result deserves commenti 


in principle, the synchronous reactance can ha 
no effect on the transient stability. 
condition 4 of section 5.2, we calculate with x 
in both axes, there can be no dependence on th: 
synchronous reactance. As was proved with referenc 
to Fig. 7, however, the stipulation of constant flu” 
linkage gives a much closer approximation to th 


{ 


true curve. With constant flux linkage the activ 
power output of the machine is an unambiguo ‘ 
function of the rotor angle 6. Fig. 9 shows this fune 
tion for a definite case in the transmission syste! 
investigated in this chapter, where the line leng 
is 300 km and the transmitted power | per uni 
The curves show the actual output of the machi 
during the various system states in terms of 6. 
curves were calculated first for x, = 0-75 and afte 


wards for x, = 2. 


A, and A, are working points of the machi 
before the short circuit occurs. The value 1 — 
is the power tending to accelerate the machine. TI 
higher the curves are situated and, in particular, tl) 
larger the included area, the better the transi 
stability of the machine, as will easily be realize 
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Fig. 8. — Curves showing the active power Py, which can be 
‘stably transmitted in the event of a three-phase short circuit with 
three-phase auto-reclosure, shown in terms of the synchronous 
reactance x, for a two-circuit transmission system 


Short circuit at the middle of the line. Line lengths 100, 200, 
300 km. 


Duration of fault =0°15 s 
Dead time =0°3s 


jThe superiority of the machine with the smaller x, 
}curves 1) stands out in Fig. 9. At this point it may 


oe underlined that the corresponding curves for the 
alculation with x’, in both axes of both machines 
ould have the same shape. 


0.3.2 The gain in stably transmissible power as a result of 
Phree-phase rapid auto-reclosure 


Fig. 10 reproduces the results of calculations with 
the same arrangement as in Fig. 9. For a line length 


pf 300 km the stability limit was computed once with 


and once without rapid auto-reclosure following inter- 
ption of the affected line, as well as with and with- 
ut voltage regulation. A striking feature is that in 
this case the gain in stability resulting from auto- 

eclosure is relatively small. Here reference may be 
n. to two publications in which analogous 
iscoveries were made ( [5], Fig. 2b on page 518; 


‘ 


0 
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Fig. 9. — Explaining the dependence of the transient stability on 
the synchronous reactances x4 and x, 
Active power output Py as a function of the rotor angle 0, 
assuming constant flux linkage. Line length 300 km. 
%q=0°65 X 0°75 p.u. 
%q=0°65 X 2:00 p.u. 


V=*%7=0-75 p.u. 
aera 200 Dy tts 


a= Normal operation 
6 = Short circuit 
c = One line disconnected 


[6], page 152). It is not surprising that the voltage 
regulation has a certain influence on the stability, 


following the remarks in 5.2. 
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Fig. 10. — Gain in stability from rapid auto-reclosure 


Line length 300 km, otherwise the conditions are the same as in 
section 5.3.1 


1 = With rapid auto-reclosure 

2 = Without rapid auto-reclosure 
a= With voltage regulator 

b = Without voltage regulator 


The remaining notation is the same as in Fig. 8. 


Conclusions 


By means of the program described information 
can be quickly obtained, within the scope of its 
capabilities, regarding a wide variety of problems 
relating to the transient stability of synchronous 
machines. This is made all the more important by 
the fact that the stability is dependent on a large 
number of factors, which make it difficult to lay 
down any general rules regarding the influence of 
different parameters. The effect of auto-reclosure, 


the synchronous reactances and the voltage regula- 


tion, as discovered by reference to various example 
cannot be transferred to quite different conditio 
without due regard to the conditions; reliable infor- 
mation can only be obtained by performing the 
respective calculation. Within the limits to which the 
program is subjected, calculation by the digital com- 
puter is clearly superior to that on the analogue 
computer, despite the fact that it takes roughly 20 
times longer to calculate one of the curves in Fig. 7 
than on the analogue computer. The deciding factor 
is the much shorter time spent in preparation. How- 
ever, if great value is placed on optimum setting of 
the regulator, the well-known advantages of the 
analogue computer become more pronounced. But | 
no particular importance was attached to this point 

. 


in the calculations described. 
W. FREY 


(KME) P. ALTHAMMER 
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